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ABSTRACT 


In the Command and Control mission, new technologies such as ‘sensor fusion’ 
are designed to help reduce operator workload and increase situational awareness. This 
thesis explored the tracking characteristics of diverse sensors and sources of data and 
their contributions to a fused tactical picture. The fundamental building blocks of any 
sensor fusion algorithm are the tracking algorithms associated with each of the sensors on 
the sensor platform. In support of this study, the MATLAB program Jusim ’ was written 
to provide acquisition managers a tool for evaluating tracking and sensor fusion 
algorithms. 

The- fusim program gives the user flexibihty in selecting: sensor platforms, up to 
four sensors associated with that platform, the target types, the problem orientation, and 
the tracking algorithms to be used with the sensors. The fusim program was used to 
compare tracking algorithms in a multiple sensor/multiple target environment. 
Specifically, the Probabilistic Data Association Filter, the Interacting Multiple Models 
Filter, the Kalman Filter and the Constant Gain Kalman Filter were evaluated against 
multiple maneuvering, non-maneuvering, and fixed targets. It is recommended that this 
study be continued to evaluate advanced tracking and data association techniques, to 
expand the program to allow attribute tracking and identification, and to study the 
Human-Machine Interface aspects of sensor fusion. 
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I. 


INTRODUCTION 


A. BACKGROUND 

The mission of airborne command and control is much more than the advertised 
“airborne early warning (AEW).” The new buzzword for the E-2C Hawkeye Command 
and Control (C^) mission is Battle Management. New technologies are being designed to 
help reduce operator workload and increase situational awareness. Ideally, these new 
capabihties will help tactical aircrews manage Battle Group assets and ultimately, the 
entire battlespace. The E-2C aircrew performs functions that cannot be rephcated by a 
machine, that is, making battlespace assessments and taking charge of situations. The 
aircrew must have a thorough understanding of the tactical picture to formulate 
appropriate responses. The industrious aircrew will endeavor to formulate the “Big 
Picture.” They mentally fuse all available information to create sufficient situational 
awareness (SA), to understand the tactical picture, and to decide on courses of action. 
However, situational awareness is easily disrupted by immediate threats, aircraft 
emergencies, or information overload. It is left up to the diligent aircrew to achieve and 
maintain SA. 

Initial detection of targets within the E-2C aircraft occurs with the radar, the lEE 
detection system, and the Passive Detection System (PDS). These systems are coupled 
with an automatic tracking system within the mission computer. Secondary sources of 
information consist of the onboard data links including, CEC, JTIDS, Eink-11, and Eink- 
4. Other data sources include SATCOM data and information received via voice 
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communications (UHF SATCOM, HF/VHF/UHF LOS). Additionally, the aircrew 
receives data such as the Air Tasking Order, Special Instructions, Airspace Restrictions, 
Rules of Engagement, Battlefield Situation Reports, and intelligence briefings prior to 
takeoff. All available information is used by E-2C aircrew to formulate battlespace 
assessments, and to create and maintain situational awareness. In a broader sense, the 
objective of any military engagement is to defeat the enemy with decisive, overwhelming 
force. Eor our military forces to take full advantage of the information and observations 
available, they must be able to interpret their battlespace and make timely and relevant 
decisions. 

Warfare technology has continued to advance, and the battlespace has become 
much more comphcated. Eurthermore, the E-2C aircrew has been increasingly tasked 
with a wider variety of more demanding missions. To keep pace with this evolutionary 
cycle of weapons advances and more complicated missions, the E-2C weapon system is 
continuing to improve. Current weapons system improvements for the E-2C include an 
upgraded Electronic Surveillance Measures (ESM) system, a new mission computer and 
displays, the Cooperative Engagement Capabihty (CEC) system, and a Satellite 
Communications (SATCOM) data system. Euture upgrades include the Radar 
Modernization Program (RMP), the UHE Electronically Scanned Array (UESA) antenna, 
the Surveillance Infrared Search and Track (SIRST) system, and the E-2C Multiple 
Source Integration system. These weapons system improvements are designed to 
improve the detection and processing of tactical information to more fuUy support the 
mission and ultimately the Battle Group Commander. 
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B. PURPOSE 


This thesis explored the tracking characteristics of different sensors and sources 
of data and their contributions to a “fused” tactical picture. This ‘fusion’ aspect implies 
that one track will represent each contact of interest in the battlespace with all available 
position, identification and intentions data displayed for each track. How can ‘sensor 
fusion’ help with formulating a tactical picture? What are the unique attributes of the 
different sensors and data sources that can be combined to produce an accurate tracking 
picture? Finally, how can the different data sources be combined to best exploit the 
attributes of identification and position information? The purpose of this thesis was to 
explore the mihtary apphcation of Sensor Fusion technology by creating a MATLAB- 
based multiple-sensor, multiple-track simulation, and exploring some of the applicable 
tracking and fusion algorithms relevant to current and future sensor fusion capabihties. 

C. TOPICS RELEVANT TO SENSOR FUSION 

I. Definitions 

First, a few definitions are necessary to ensure that no terminology problems are 
encountered throughout the remainder of this paper. Currently, sensor fusion takes many 
forms, whether from industry, government services, or military apphcations. The 
primary focus of this hst is the mihtary application of sensor fusion. 

Battlespace 

The term “battlespace” includes all of the area (land, sea, air, and space) of interest, ah of 
the sensors and information sources that can contribute to a coherent tactical picture, ah 
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of the forces (friend, enemy, and neutral), weapon systems, and aU communications 
systems in the area of interest. 

Kinematic Data 

Target kinematics are usually represented by the position, course, speed, altitude, or 
accelerations and turn/climb rates [Ref. 1]. 

Attribute Data 

The target attributes are the features of the target such as returned signal strength, IR or 
radar signature, or target classification and ID information [Ref. 1]. 

State 

The target state [Ref. 2] is a determination of target position, course, and speed, 
depending on the current sensor measurement and the previous state estimate. The state 
vector usually consists of the target position coordinates and the corresponding velocity 
components. The target state can also hold estimates of turn rates or other kinematic 
data. 

Local Active Track 

A local active track (LAT) refers to a track held by a sensor that is organic to the sensor 
platform. A track [Ref. 1] is the result of a state prediction from the tracking algorithm, 
and is based on the measurements from the sensor. 

Similar Source Integration 

Similar source integration (SSI) is the association and correlation of tracks from the same 
kind of sensors such as radar or IFF detection systems [Ref. 4]. 
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Dissimilar Source Integration 


Dissimilar source integration (DSI) is the assoeiation and eorrelation of tracks from 
multiple types of sensors such as an IFF-to-radar correlation or a LAT-LAT/Data Link- 
LAT correlation. A DSI must include the sensor registrations for time eorreetions and 
coordinate transformations [Ref. 4]. 

Multiple Source Integration 

The multiple source integration (MSI) process performs the data association and 
correlation of tracks from multiple types of sensors. The output of the MSI proeess is a 
database of aU possible track-to-track associations [Ref. 4]. 

Data Fusion 

The most recent definition of “Data Fusion” is from the Joint Directors of Laboratories 
(JDL) paper of 1998. According to the JDL [Ref. 3], “Data Fusion” is the process of 
combining data to refine state estimates and predictions. 

Sensor Fusion 

Strictly speaking, the term sensor fusion is the combination of estimates from multiple 
sensors for the formulation of multiple tracks [Ref. 1]. The term “Sensor Fusion” is 
thrown around loosely and is sometimes taken to mean the same as the ‘data fusion’ term. 
While data fusion is truly a more encompassing term, the ‘sensor fusion’ term imphes a 
distinct relationship between the measurements/observations of the battlespaee sensors 
and the fusion of the resulting information. Therefore, ‘sensor fusion’ will be used 
throughout the remainder of this paper. 
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2. Combat Identification 


Identification is the key element in the decision making process, i.e., the 
identification of a contact leads to a decision. Given the identification of a contact, the 
decision-maker can decide whether to investigate, monitor, ignore, intercept, or shoot the 
contact. For example, if a contact is flying in a zone where no friendly aircraft are known 
to be flying (with 100% certainty), the specific identity or intentions of the target may be 
irrelevant and stiU satisfy the conditions of combat identification (CID). Deciding which 
contacts to process (and when) is a function of the threat (situation) assessment and the 
assets available to meet the threat (impact assessment). Situation and impact assessment 
wiU be discussed in the next chapter. 

The military apphcation of sensor fusion is primarily the fusion of all organic 
sensor track data, off-board data hnk information, and archived information to form a 
single fused track for each contact of interest in the area of interest. In the CID process, 
three vital pieces of information about a contact are its position, its identification, and its 
intentions. Depending on the source of information, these attributes can be closely 
related, such as a radar track with IFF information fused by a tracking system, or as 
diverse as an ESM hue of bearing that cuts through a postulated ground emitter site. The 
diversity of this data is precisely the purpose for studying the techniques of sensor fusion. 
The use of ID information leads to the ultimate goal of combat identification (CID) for 
weapons system engagement of enemy targets. 
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3. The Operator’s Point of View 

The goal of sensor fusion is to ease the workload and improve the situational 
awareness of the crew. Sensor Fusion systems are meant to aid the crew in the decision 
making process. In other words, let the machines do the number crunching, and let the 
humans make the decisions. One of the most difficult questions to answer is: What does 
the operator want? Unfortunately, an old adage stiU apphes: If 20 different operators are 
asked for their views on how information should be displayed, the result will be 20 
different answers. The operator is the primary interface between the display and the 
command and control decision tree. In the case of the E-2C, each operator can act 
independently. 

With respect to sensor fusion, the individual display of information must address 
the diversity of missions within the E-2C or any other command and control platform. 
How the information is displayed is vitally important. The display of information can be 
prioritized depending on the tactical situation and the tactical impact of the target. The 
operator that is controUing a surface search coordination (SSC) mission does not have the 
same information and threat display priorities as the operator who is providing tactical 
control for an overland strike mission, even though both events may be happening 
simultaneously. Instead, the operators need maximum flexibility for prioritizing their 
individual display requirements for information, threats, timehnes, etc. Eor the airborne 
battlefield command and control mission, displaying all mission-apphcable information 
would be senseless for the user and impossible to fully interpret. Instead, the display of 
information should reflect a relative importance to the mission at hand, the desires of the 
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operator, and the needs of the battlespace commander. Ideally, only those threats that are 
most apphcable to the operator’s mission would be displayed. 

D. CHAPTER OUTLINE 

Chapter II is a discussion of Sensor Fusion and why it is important to the mihtary. 
Chapter III is a description of modern tracking algorithms. Chapter IV is a detailed 
description of the MATLAB based fusion simulator program, including the evaluation 
parameters of the test sets. Chapter V describes the results of testing the tracking 
algorithms for stabihty, and presents the results of running several different test cases. 
Conclusions and recommendations are included in Chapter VI. 
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II. OVERVIEW OE SENSOR EUSION 


A. WHY SENSOR FUSION? 

Simply put, tactical aircrews have been doing sensor fusion all along by achieving 
situational awareness with respect to the tactical picture. This thesis was written with the 
platform and the operator as primary concerns. The importance of the operator and the 
platform cannot be overstated: every platform has unique sensor capabihties and 
attributes, communications systems, and display capabihties. And every operator has 
unique display requirements, depending on the platform, the sensors and data sources, 
and the mission being conducted. The sensor/comm suite and display capabihties of the 
AEGIS cruiser and the E-2C are vastly different, yet the C goal is still the same. The 
operator must be able to interpret and respond to the tactical situation depending on 
his/her mission goals. Accurate tracking and identification are paramount to this process. 
Primarily, the operator would hke to know the identification of a contact so the tactical 
impact can be determined. The operator must consider the impact to his platform, the 
platforms under his control, the carrier, and ultimately, the entire battle group. This is 
how Combat Identification becomes an important by-product of the sensor fusion ID 
process. If engaging a contact is within the ROE, the operator must be able to ahocate 
assets quickly and decisively without distractions from information overload, false tracks, 
or dual tracks. 

Instead of cluttering up the scope with thousands of reports, the information 
would ideally be sorted by target type, registered in time to allow a common reference 
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time, and fused with other reports and local active tracks. This data association provides 
an identification of a track, or increases the confidence of the position report. Even if 
there is no association with a local track, such as in the case of a known SAM site, the 
information must be displayed in a timely and relevant manner. Additionally, the 
information could be used to cue the operator to look for contacts in an area where no 
active track exists, or to cue the radars of other platforms to search the area. Contact 
reports that are considered an immediate threat should pop up and alert the operator so 
that the operator can respond appropriately. The display of threat information should 
coincide with the operator’s display desires and the tactical situation. For example, 
certain high-priority threats can be displayed anytime, anywhere on the scope. Other 
threats/tracks may only be displayed in an area of interest, depending on the mission or 
the needs of the operator. 


B. THE LEVELS OF SENSOR FUSION 

The Joint Directors of Laboratories (JDL) data fusion process model is the most 
commonly used model to communicate ideas about algorithms, systems, and research 
(Figure 1). In its everyday use, the term “data fusion” generally refers to that of Level 1 
fusion. This is not the most accurate interpretation of the term. In Level 1 data fusion, 
all sources of information, including track data and attribute data, are combined to form a 
composite tracking picture. In order to more fuUy understand the big picture on sensor 
fusion, all 5 levels are discussed below, starting with the Data Fusion Functional Model 
[Ref. 5]: 
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Figure 1. The Data Fusion Functional Model*. 
* Reproduced From [Ref. 5] 


1. Level 0 - Sub-Object Data Assessment: 

Level 0 is considered to be pre-object processing. By definition [Ref. 3], Level 0 
fusion is, “estimation and prediction of signal/object observable states on the basis of 
pixel/signal level data association and characterization.” In other words, what kind of 
targets or signals are expected in the area of interest? From this early assessment, we can 
estimate and predict the signal or object observable states. Types of data would include 
signal detected on the basis of integration of a time-series of data or feature extraction 
from a region in imagery. 


2. Level 1 - Object Assessment: 

Level 1 fusion is the level that is normally associated with the term “sensor 
fusion,” and is the primary focus of this Thesis. At this level, the raw measurements are 
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used to estimate the current states of each entity. AU organic sensor information is 
combined with off-board information including data hnks, SATCOM data, Order of 
Battle databases, ATO/SPINS/ROE, and other database information. During “object 
assessment,” estimation and prediction of object states occurs based on the measurements 
and/or observations. The kinematics and attributes of the object answer the questions of 
where and who the object is. According to the JDL [Ref. 3], Level 1 fusion is, 
“estimation and prediction of entity states on the basis of observation-to-track 
association, continuous state estimation (kinematics), and discrete state estimation (target 
type and ID).” Mainly, Level 1 fusion provides the composite tracking picture. 

3. Level 2 - Situation Assessment: 

In Level 2, a thorough threat assessment is conducted based on the current 
inteUigence reports and other measurements of troop/hardware movements, 
communications, etc. By JDL definition [Ref. 3], Level 2 fusion is, “estimation and 
prediction of relations among entities, to include force structure and cross force relations, 
communications and perceptual influences, physical context, etc.” Level 2 fusion 
involves associating the hypothesized entities or tracks into aggregations. According to 
the JDL [Ref. 3], the aggregate state can be represented by a network of 
interconnectivity and relations between the elements. The relations considered include 
physical, organizational, informational, and perceptual according to the mission of the 
entity. The term ‘situation’ is used to describe an aggregate object of estimation. 

There are many different approaches to fusing the information at Levels 2 and 3. 
Lor example, Bayes theorem and Bayesian behef networks, neural networks, fuzzy logic. 
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genetic algorithms, and Hidden Markov Models (HMM) [Ref. 5] can be used to model 
the states and estimate situations. 

4. Level 3 - Impact Assessment: 

The JDL [Ref. 3] definition of Level 3 fusion is, “estimation and prediction of 
effects on situations of planned or estimated/predicted actions by the participants to 
include interactions between action plans or multiple players (e.g., assessing 
susceptibilities and vulnerabihties to estimated/predicted threat actions given one’s own 
planned actions).” The Impact Assessment is also known as Threat Refinement [Ref. 6]. 
Threat Refinement is used to estimate enemy capabihties, identify threat opportunities, 
estimate enemy intent, and determine levels of danger to the friendly forces. In doing so. 
Level 3 fusion focuses on estimating the l ik elihood of hostile actions, and determining 
the projected outcomes if hostile actions do occur. Part of Level 3 fusion [Ref. 5] is 
threat prediction, which attempts to answer the following question: Based on enemy 
capabihties (firepower and preparedness), what enemy actions pose a threat and to what 
extent could the enemy damage friendly forces? One key element to predicting the threat 
is determining the possible enemy courses of action (COA). These COA’s could be 
based on lessons learned form previous actions, or hypothesized based on terrain and 
natural barriers such as cliffs, lakes, rivers, and oceans. To quantify the COA’s, the 
likelihood of each possible CO A is estimated based on the observations. 

During the generation of COA’s, the information from the lower level fusion 
engines would be necessary for accurate and sensible predictions. For each course of 
action, the impact to blue forces must be assessed to determine the probable outcome. 
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From this evaluation, the most threatening COA’s can be identified for the battlefield 
commander. For example, the Enemy Order of Battle would include all of the enemy 
forces in the battlespace. Based on the units involved, the hardware, the terrain, and 
many other observations about the enemy, the candidate COA’s can be determined. 
Depending on their movements (the measurements or observations), the li kelihood that 
the movements are associated with a particular COA can be determined. This l ik elihood 
can be determined for each candidate COA, and the unlikely COA’s can be discarded. At 
this point, the possible outcomes of particular enemy actions can be determined. Are the 
enemy forces assembhng, attacking, defending, or deploying? Or, is a particular unit 
attacking while another is staying back to defend against a counter attack? Next, a 
determination can be made about the likely outcome of an engagement for each possible 
COA. From this determination, the threat level to friendly forces can be predicted using 
a simple decision tree. Finally, an overall hkely threat to friendly forces from all possible 
enemy courses of action can be determined. Should an engagement occur, what is the 
likely outcome? 

5. Level 4 - Process Refinement: 

First, a definition is necessary to complete this section. According to JDL [Ref. 
3], Level 4 fusion is, “adaptive data acquisition and processing to support mission 
objectives.” In Level 4 fusion, resources are assigned tasks to support the formal 
relationship between estimation and control, and association and planning. Simply put. 
Level 4 fusion is resource management to support the particular operator needs. Lusion 
across the levels takes place to arrive at the most probable enemy course of action based 
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on the measurements, the enemy’s capabihty, and the possible courses of action [Ref. 7]. 
The operator requirements drive the optimization process of Level 4 fusion. Particular 
information needs can be used to cue sensors or to optimize the fusion process at each of 
the other levels. 
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III. SENSOR FUSION TECHNIQUES AND PRACTICES 


A. TRACKING FUNDAMENTALS 

1. General 

The fundamental building block of any sensor fusion algorithm is the tracking 
algorithm associated with each of the sensors on the sensor platform. Issues relevant to 
tracking are: track initiation, track maintenance, drop track criteria, and handling of false 
tracks. Additionally, how the raw data is combined within the processors and mission 
computer can have an impact on how the data is combined at the data fusion level. 
Primarily, this discussion will focus on the algorithms required for track ‘maintenance’ 
vice track initiation or drop track criteria. 

2. Sensors and Sources 

In an ideal situation, every threat and possible outcome of an engagement could 
be considered prior to firing a first shot. Naturally, the quahty of fusion information 
given to decision makers is only as good as the sensors and sources of the information 
that are used as inputs. Sensor types are as diverse as the platforms that carry them. 
Radars can vary by mission, type, low/high frequency, scan types, pulse repetition 
frequencies, and by many other differences. Electronic Surveillance Measures systems 
vary in their abilities to distinguish signal parameters, measure angles of arrival, or 
determine signal identification. For the purposes of Sensor Fusion, sensors and sources 
provide the raw and preprocessed data to fusion systems for composite processing. The 
following definitions are included in Waltz and Fhnas [Ref. 2]. 
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a. Sensors are devices that detect and/or measure physical phenomena. 
Remote Sensing describes a sensor that measures a distinguishing 
phenomenon that was transmitted through a medium, such as with an 
ESM system. 

b. Sources describe the variety of originators of data, such as 
observations, intercepted communications or other data, a priori 
information such as map data (terrain, roads, cities, lakes, rivers, etc.), 
sea lanes, flight routes, etc., and other archived data such as the 
OOB/EOB, the ATO/SPINS/ROE, and intelligence data. 

c. Einks are the communications and connectivity from the sensors and 
sources to the data fusion processing nodes. 

3. Detection and Tracking 

Sensors detect the presence of signals and estimate the parameters of those 
signals. The signal processing of the sensor data results in a measurement for that current 
sample time [Ref. 1]. Tracking is the result of processing the measurements for the 
purpose of estimating the current state of the target. The state estimate can consist of the 
target kinematics and attributes. According to Bar-Shalom [Ref. 1], measurements are 
the observations of the target, and are usually corrupted by noise from the processing 
sensor or the signal transmission medium. The characteristics of the measurement are 
dependent on the sensor type making the measurement. Eor example, a 3-D radar would 
directly measure the range, azimuth, and elevation to the target (relative to the sensor). A 
2-D radar would measure the range and azimuth to the target. A passive system such as 
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an ESM system would measure the direction of arrival (DOA), but could also determine 


many of the signal characteristics (such as signal strength, frequency, and pulse repetition 
frequency). The noise in the measurements is the result of uncertainty in the process, 
such as detection false alarms, clutter, other targets, and deception/countermeasures [Ref. 
1]. For the purposes of this thesis, track initiation and drop criteria are not addressed. 


4. Data Association 

Any discussion of tracking algorithms would be incomplete without a description 
of data association methods. While the term ‘tracking’ is used for describing the state 
estimation process, ‘association’ is used to indicate the process of matching the data or 
knowledge/information to the track [Ref. 1]. The data association process answers the 
question, “which set of data belong to which track?” According to Waltz and Lhnus 
[Ref. 2], data association is, 

the process of relating individual sensor measurements (data) to other 
measurements to determine if they have a common source (e.g., target or event). 
Although the measurements may be referenced to different coordinate systems 
with different spatial accuracy’s or resolutions, the association process must relate 
each measurement to a number of possible sets of data, each representing a 
hypothesis to explain the source of the measurement: (1) The false alarm set, 
indicating that the measurement is unreal and to be ignored, (2) The new target 
set, indicating that the measurement is real and relates to a target for which there 
are no previous measurements, (3) An existing set of previous measurements 
related to a single target. A set exists for each previously detected target. 


By Bar-Shalom [Ref. 1], the three categories of data association are measurement-to- 
measurement association, measurement-to-track association, and track-to-track 
association. Measurement-to-measurement association is used for track initiation. 
Measurement-to-track association is used for maintaining the track, such as predicting the 
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current state of the target based on the hkelihood that the measurement is actually from 
the target. Track-to-track association is used in a multiple sensor environment, when two 
or more sensor systems are reporting weU-estabhshed tracks. When measurements (or 
track data) are available from more than one sensor, a process called ‘data registration’ is 
required to ahgn the samphng times of the sensor data before track-to-track association 
can be used to combine the data [Ref. 1]. 

The methods of data association require a li kelihood or probability measure to 
evaluate alternatives [Ref. 1]. For example, a measurement could be ‘associated’ with a 
current track, another track within the same detection region (gate), a clutter point, or 
some other anomalous situation. Since it cannot be known with absolute certainty what 
the measurement is actually associated with, the previous hypotheses or guesses are 
typically evaluated based on the lik elihood that they are the correct guess. Typically, 
target gates are used to eliminate measurements from consideration in the development of 
hypotheses. Measurements that fall outside of the gated region are not used in the state 
estimation of the target. Target ‘gates’ are discussed in the next section. Only the 
measurements that fall within the gate are considered in estimating the state of the target, 
depending on the algorithm used. Each measurement can be assigned a lik elihood or a 
probabihty that it is the measurement that is most correct for updating the state. 

Numerous statistical approaches exist for assigning li kelihoods or probabihties to 
the data association process. Several of the methods for data association include the 
Maximum Likelihood Estimation, Bayesian Statistics, Evidential Reasoning, Dempster's 
Rule, and Dempster-Shafer Reasoning. 
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5. Gating 

Gating is used to narrow the search around a predicted target position for the next 
update or measurement. A three dimensional region or volume of probability surrounds 
the predicted position of the target. According to Blackman and Popoh [Ref. 8], gating is 
used to eliminate unhkely observation-to-track pairings. If an observation falls within the 
gated region, the track can be updated with that observation. Measurements or 
observations that fall outside the gate are not considered in updating the target state at the 
time of the measurement. Which returns within the gated region will be used to update 
the track depends on the data association algorithm used. Data association methods will 
be discussed in the next section. A gated region can take many different forms including 
rectangular, eUipsoidal, spherical, etc. 

Maneuver gating is also an important consideration, depending on the tracking 
algorithm used. According to Blackman [Ref. 8], a maneuver gate closely models the 
most severe potential maneuver that a target can perform. Accordingly, the region can be 
expanded or stretched in any direction depending on the anticipated maneuver. 

Three simple approaches to data association within a tracking gate follow: 

a. Nearest Neighbor 

The Nearest Neighbor (NN) approach is the simplest measurement-to- 
track association method. With NN, the measurement that is closest to the predicted 
position of the track is used to update the track state, provided that the measurement is 
within the tracking gate [Ref. 1]. 
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b. Global Nearest Neighbor (GNN) 

The global nearest neighbor (GNN) method is a very simple and widely 
apphed method for data association [Ref. 8]. Under the GNN method, a hypothesis for 
aU the possible associations is made for updating the target state. Only the most lik ely 
association is used to update the state and the other possibihties are discarded. Problems 
occur when more than one measurement occurs in the tracking gate and the formation of 
a solution occurs with an association assignment matrix. 

c. All Neighbors Method 

The All Neighbors approach is used in the Probabihstic Data Association Filter, 
which will be discussed in the section on tracking algorithms. 

B. TRACKING ALGORITHMS 

1. Batch Processing 

Numerous techniques have been used for this simplest case of tracking a single 
target that is not maneuvering. The approach is to collect a number of hits on the target 
and batch process the data to produce a track [Ref. 9]. With batch processing, the more 
hits collected, the better the solution. However, the computational requirements are more 
demanding as the number of hits increases, making batch processing an impractical 
solution for most surveillance systems. Every time a new measurement is collected, all 
the previous measurements are used to calculate the new state, making batch processing 
very cumbersome. 
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2. Prediction-Correction Methods 


One of the advantages of prediction-correction, or recursive, methods is that the 
state update depends only on the previous state and the current measurement. According 
to Hutchins [Ref. 9], the simplest trackers are only designed for straight-hne motion 
(SLM), yet target accelerations can be handled by increasing the gain of the filter or by 
increasing the noise term of the state equations. The more comphcated trackers are 
designed to compensate for turning motion, or to handle multiple hypotheses of target 
motion, or to sort through clutter. Tracker types can be combined to allow multiple 
sensor, multiple target tracking, and data association. 

a. Alpha-Beta Tracker (7-7) 

The Alpha-Beta tracker is the simplest constant-gain tracking algorithm. 
The tracker has poor performance, but requires a very low computational load. In this 
case, a constant gain matrix is set up for the target position update equation. The alpha- 
beta tracker is used in tracking systems where position measurement updates are 
available and the state vector consists of positions and velocities [Ref. 9]. The value of 
the gain is preset for handling straight-hne motion or turning motion. When the gain is 
set to compensate for turning motion in a target, the straight-hne performance wih suffer 
somewhat. This is also true for many of the trackers in the discussions that follow. An 
extension of the alpha-beta algorithm is the alpha-beta-gamma tracker, which includes 
accelerations in one state vector [Ref. 8]. 
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b. Constant Gain Kalman Filter (CGKF) 


The Constant Gain Kalman Filter (CGKF) is a simplified case of the 
Kalman Filter [Ref. 9]. Instead of updating the covariance matrix at every measurement 
update, the covariance is taken to be constant. Using the assumption that the covariance 
matrix wiU approach a steady state value over time, the Algebraic Riccati Equation 
associated with the hnear, time-invariant, discrete time system can be solved numerically. 
For this case, the MATLAB function dlqe (discrete-time linear quadratic estimate) can be 
used to determine the values of constant covariance and constant Kalman gain. Then the 
Kalman update equations depend only on the previous state, the measurement, and the 
constant gain. Although the CGKF is not an optimum solution, it is not heavily burdened 
computationally. 


c. Kalman Filter (KF) 

The Kalman Filter (KF) is the basis upon which many of the more 
advanced algorithms are designed. The Kalman filter is an optimum solution to the 
sequential least squares problem, meaning that the least square error is minimized. It is a 
sequential algorithm because it only depends on the most current measurement and state 
estimate, and the associated measurement and state prediction covariance matrices [Ref. 
9]. The Kalman filter is not very computationally demanding, but the filter is not 
designed to handle maneuvering targets, clutter, or multiple targets. The filter can be 
adapted to handle maneuvering targets, but the solution is no longer optimum. Examples 
of this phenomenon wiU be shown in the Results section. 
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d. Extended Kalman Filter (EKF) 


The Extended Kalman Filter (EKF) is used in cases where the mapping of 
coordinates is non-linear. The EKF is most apphcable in cases where the measurement 
process is non-linear, or the target dynamics are non-hnear [Ref. 8]. According to Bar- 
Shalom [Ref. 1], the nonhnear transformation may introduce bias in the solution, the 
covariance calculation is not necessarily accurate, and the EKF can diverge if the initial 
conditions are inaccurate. The EKF was not examined thoroughly in this thesis. 

e. Interacting Multiple Models (IMM) 

The Interacting Multiple Models tracker is used to predict the current state 
of the target using two or more different models. For example, if the target is expected to 
be a maneuvering target, the models used could be a straight-line motion (SFM) model, a 
left turn model, and a right turn model. Other models used could be a variety of turn rate 
models or climbing/descending models. The number of models used is dependent on the 
apphcation. In this thesis, the IMM used is a 2-model IMM, where the two models differ 
only by the noise term (one for SFM, and one for turning motion). 

For the IMM estimator, multiple state equations are used to describe each 
of the different modes of operation [Ref. 9]. A Markov transition matrix is used to 
specify the probability that the target is in one of the modes of operation. Usually, these 
values are chosen heuristicaUy. In the case of the 2-model IMM used in this thesis, the 
chosen probabilities were (1) A 10% probabihty that the target would turn if in SFM at 
the measurement time, and (2) A 33% probabihty that the target would return to SFM if 
in a turn at the time of the measurement. Similar to the “soft-switching” discussed in an 
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Air Force Research Lab report [Ref. 10], the model probabihties are updated at each new 
measurement, and the resulting weighting factors are used in calculating the state. In 
other words, no gating decision is required for the tracker to function properly. See 
Chapter IV for a more thorough description of the algorithm. 

/. Multiple Hypothesis Tracker (MHT) 

The MHT is a very complex, yet flexible approach to solving the multiple 
target data association problem [Ref. 9]. Even in a very simple case where one 
measurement can be associated with one track, numerous possible hypotheses exist as to 
the precise nature of the association: (1) the measurement may be associated with the 
currently held track, (2) the measurement may be associated with a new track, (3) the 
measurement may be associated with no track. With multiple tracks and multiple 
measurements, the number of hypotheses can grow dramatically. Each hypothesis 
receives a score, which is dependent on the probabihty (hkehhood function) that the 
hypothesis is correct [Ref. 9]. AU the hypotheses can be maintained from scan to scan, or 
the hypotheses can be pruned back to the most likely one or two track solutions. Eor this 
reason, an MHT is a very computationally demanding algorithm and it requires a great 
deal of memory. 


g. Probabilistic Data Association Filter (PDAF) 

The Probabilistic Data Association Eilter (PDAE) is used to compensate 
for clutter that is received as possible vahd target returns within a gated region around the 
predicted target position. The PDAE is a suboptimal Bayesian algorithm that associates 


26 



probabilistically (using the likelihood function) aU the vahdated measurements to the 
target of interest [Ref. 1]. For this reason, the PDAF is also known as the “all-neighbors” 
data association approach. The validated measurements, also known as the neighbors, 
are combined using a weighted li kelihood function in the algorithm to account for 
measurement uncertainty. The state estimate is updated with aU the validated 
measurements weighted by their lik elihoods of having originated from the target, i.e. the 
combined innovation. To account for the measurement uncertainty, an additional term is 
added to the covariance update equation. As opposed to the nearest neighbor approach to 
data association, the PDAF is an ‘all neighbors’ data association approach. 

h. Joint Probabilistic Data Association Filter (JPDAF) 

The Joint Probabilistic Data Association Filter (JPDAF) is an extension of 
the PDAF to include multiple targets in clutter [Ref. 1]. The JPDAF allows overlapping 
vahdation regions, meaning that more than one target may be present within the gate for 
each update. This effect causes a persistent interference over several samphng times. 
One of the simplifying assumptions of the JPDA approach is that the number of targets is 
known. Also, the false measurements are uniformly distributed across the vahdation 
region, meaning the extra computation time is expended evaluating all measurements as 
if they were vahd. 
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c. 


ADVANCED TRACKING ALGORITHMS 


1. General 

In order to handle multiple sensors and multiple targets that may be maneuvering, 
additional methods are necessary. These methods improve on the basic characteristics of 
the trackers discussed in the previous section. 

2. Multiple Targets 

The multiple target tracking problem is discussed extensively in Bar-Shalom [Ref. 
1]. If the number of targets is know, the JPDAF is a weU estabhshed algorithm for 
tracking the targets. Since the JPDAF is a Bayesian approach, all the possible targets that 
can be updated by a measurement are considered simultaneously. Additionally, aU the 
vahdated measurements that faU within the tracking gate are considered in updating the 
track. However, if the targets are maneuvering, or the number of targets are unknown, 
some other implementation of the filter is required. 

Advanced tracking algorithms are designed to combine the best attributes (or the 
most practical attributes) of the more basic tracking and data association techniques. For 
example, combining the IMM with the JPDA allows multiple model tracking with the ‘all 
neighbors’ approach to data association. In fact, many of the algorithms can be 
simplified to allow a sub-optimum solution while preserving the practical aspects of the 
algorithm. For example, the JPDA could be greatly simphfied by forcing nearest 
neighbor data association (NNJPDA). Combining this new filter with an IMM would 
result in a filter that could track maneuvering targets with nearest neighbor data 
association (IMMw/NNJPDAF). 
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3. Multiple Sensors 

In multiple sensor tracking scenarios, the sensors can be of the same type in 
different locations, or of different types in the same or different locations. A primary 
consideration in the multiple sensor problem is the time at which the measurements are 
collected. If the sensors are perfectly synchronized in time, measurement-to- 
measurement associations can occur, followed by centrally located association and 
tracking [Ref. 1]. A more common approach occurs where each sensor performs the 
measurement-to-track association and tracking, and then track-to-track association and 
fusion occurs between the sensors. 

When performing track-to-track association and fusion, a test must be performed 
to determine if two tracks belong to the same actual target. This is done by comparing 
the state estimates of the two tracks [Ref. 1]. To test the hypothesis that two tracks 
belong to the same target, the difference between the state estimates is compared to a 
threshold value. If the test passes, the fusion of the estimates can be carried out with a 
simple equation that treats the error covariances of the two sensors independently. The 
resulting state estimate is usually more accurate than the single sensor case. 

D. ATTRIBUTE TRACKING/FUSION 

According to Drummond [Ref. 11], the use of the term ‘attributes’ to describe all 
the characteristics of a target is too broad due to the way that the data is treated for 
tracking purposes. By Drummond’s definition, features, attributes, and categorical 
features, are more descriptive terms. A feature is a characteristic obtained from a 
continuous sample space to include examples like target size, radar cross section, or other 
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signature data. An attribute is based on information obtained in a discrete sample space. 
Examples of attributes include target type, radar type, number of engines, or the IFF 
information. The distinction is made between features and attributes because of the way 
the information is processed. When Bayesian techniques are used, features are processed 
according to their probabihty density, and attributes are processed based on their 
associated discrete probabihties. Categorical features are another type of data that 
include information that is already known, such as the wingspan of potential target types, 
or expected IFF information. This information can be measured directly, obtained over 
time, or combined in a way that allows a comparison to known characteristics of targets, 
thereby allowing classification into a finite set of categories. 

The feature and attribute data is most commonly used for target classification and 
identification, or for Combat Identification (CID). This is an important concept in Fevel I 
Sensor Fusion, the detection and estimation of track attributes. However, a distinction is 
made between these inherent features or qualities of a target and the behavior of the 
target [Ref. 8]. The perceived target behaviors or intentions are the result of higher level 
sensor fusion, namely the situation assessment (Fevel 2). According to Blackman [Ref. 
8] and Drummond [Ref. 11], the features and attributes of a target can be detected, passed 
through a thresholding process, estimated, and tracked, similar to tracking based on target 
kinematics. 
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E. SENSOR FUSION 


The sensor fusion problem is a culmination of all the aforementioned tasks 
associated with sensor/multiple sensor tracking, attribute tracking, and track-to-track data 
association. One of the key features of a sensor or data fusion system or algorithm is the 
necessity to perform data registration. Also, in order to share the information generated 
in a multiple sensor environment, a robust communications system must be in place to 
allow timely and accurate sharing of track information. 

The sensor registration problem necessitates a correction in track databases due to 
systems that are not synchronized in time, and corrections due to common reference 
frame differences (such as gridlock problems). Also, according to Blackman [Ref. 8], 
many other sources of registration error can contribute to misahgnment in track 
information. 

1. Coordinate system errors: misahgnment of coordinate axes in the 
measurement systems. 

2. Bias error: errors due to range and bearing bias error. 

3. Location error: errors due to variation in navigation solutions. 

4. Other error sources: radar refraction/ducting errors, bias errors from 
polar-to-Cartesian coordinates. 

FinaUy, the architecture of sensor fusion systems are addressed in numerous 
publications including [Ref. 12], and are not thoroughly discussed in this paper. Whether 
a fusion system wiU process data in a centrally located database management system or 
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be distributed to the sensor platforms is an important consideration. In the case of a 
distributed system, the platform that has the most correct solution about the position, 
identification, and intentions of a contact may not be readily known. Indeed, the 
confidence in the data needs to be addressed for both types of system architectures. 

The information received may be vital to the mission. A Battle Group 
Commander may want to maximize the situational awareness of his platforms via a 
common operational picture. An aircraft responding to a time-critical strike request will 
need max im um information about target location and movement, target type or 
identification, and possible threats to his platform. 
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IV. FUSION SIMULATOR 


A. FUSIM OBJECTIVE 

The objective of the fusim program is to provide acquisition managers a tool for 
evaluating tracking and sensor fusion algorithms. Also, this first edition of the program 
provided a basis of study for this Thesis. The fusion simulator was written in MATLAB 
6.0 for the following reasons: (1) The MATLAB software is widely accessible and runs 
on any personal computer. (2) MATLAB is compatible with other languages including 
C++ and Java. (3) MATLAB is highly flexible and optimized for vector/matrix 
operations. (4) MATLAB can be programmed to utihze multi-dimensional arrays. (5) 
MATLAB can be used in an object-oriented programming environment. 

The fusion simulator (FUSIM) was written from the perspective of the sensor 
platform and the operator. Any type of platform can be selected for the sensor platform, 
either airborne or surface, maneuvering or non-maneuvering, or a fixed site. FUSIM is 
very flexible for selection of targets. The user can select any number of high-speed/low- 
speed/fixed maneuvering or non-maneuvering targets for maximum flexibility in a 2- 
dimensional simulation environment. Also, the user can select up to 4 different sensors 
for the sensor platform. For output, the user can select a real-time output plot (PPI) or the 
PPI and the sensor/target/tracker Least Square Error plots. 
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B. FUSIM DESCRIPTION 


The fusim program was written to give a user flexibility in selecting sensor 
platforms, up to four sensors associated with that platform, the target types, the problem 
orientation, and the tracking algorithms to be used with the sensors. Table 1 is a hst of 
potential/M5/m configurations that a user could select. 


Sensor 

Platforms 

Airborne 

Maneuvering 

Surface 

Moving - Maneuvering 

Non-Maneuvering 

Moving - Non-Maneuvering 

Fixed (like an aerostat radar) 

Fixed Site 

Potential 

Sensors 

Airborne 

Surveillance 

Radar 

Airborne 

Fighter 

Radar 

Non- 

Cooperative 
Target Radar 

IFF Interrogate 

Airborne ESM 
error model 

Infrared Search 
and Track error 
Model 

Surface 

Surveillance 

Radar 

Surface 

Fire 

Control 

Radar 

Eastern IFF 

IFF Reply 

Surface ESM 
error model 

Ownship 

Navigation 

model 

Target 

Types 

Airborne 

Maneuvering 

Surface 

Moving - Maneuvering 

Non-Maneuvering 

Moving - Non-Maneuvering 

Fixed 

Fixed Site 

Tracking 

Algorithms 

Probabilistic Data 
Association Filter 
(PDAF) 

Interacting Multiple 
Models (2-model) Filter 
(IMM) 

Kalman Filter (KF) 

Constant Gain Kalman 
Filter (CGKF) 


Table 1. List of Potential Configurations for the. fusim Program. 


Within the fusim program, the FUSIM file is the main control file for the entire 
simulation. FUSIM calls the various associated functions to set up the needed matrices 
and variables, creates the truth data for the targets, adds the measurement error, sends the 
noisy measurements to the trackers, and then stores and plots the target states and 
associated least square errors. The FUSIM function is called from the MATLAB 
command line. A complete description of the fusim program and the related functions is 
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included in Appendix A. The MATLAB code for aU the functions is included in 


Appendix B. 


C. THE MATHEMATICS OF FUSIM 

The fusim program runs on a few simple concepts of modehng and simulation. 
For example, the program utihzes synchronous timing for the sensors associated with the 
sensor platform, making the data association problem a httle easier to deal with. This 
section deals with the mathematics of the fusim program, describing the construction of 
the various matrices and vectors, and describing their interactions at specific points in the 
program. 

1. Startup 

Fusim begins with the START function for initiahzing a variety of variables 
including the Extraction Matrix (H) and the Discrete Time State Equation Transition 
Matrix (E). Also, aU the target matrices are initialized and the simulation timing is 
determined. 

Extraction Matrix (H) 

The Extraction Matrix simply extracts the x,y components from the 4x1 state vector. As 

used in this program, H is a 2x4. 

?x? 1x1 positionl 

7 7 7 7 

?1 0 0 0? AvxA Ax ? velocity A 

HI ^ 7 ; state 1 I l = l "I 

7OOIO7 1 yl lyl positon 1 

777 7 

7 vy 7 )yl velocity) 
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Total Simulation Time (t) 


The total simulation time in minutes is input by the user via the START.MAT file. After 

reading in the value, the simulation time is converted to seconds. 

simstart ? %imulationtime,nloops,deltatime,DLRP,N!A, filterselectionll 4? 

t ? simulationtime * 60 

Simulation Time (simt) 

The cumulative simulation time is stored in the (simt) variable. The value of (simt) starts 
at zero and runs up to the total number of seconds in the simulation. 

Number of Simulation Loops (nloops) 

The value of (nloops) is usually 1 for a simulation with a real-time output. For multiple 
runs, a value of 20 provides a good indication of performance. For sufficient averaging 
of the Monte Carlo runs, a value of 100 or 200 is required. However, the computation 
time is dramatically increased. 

Samphng Time (delta) 

The sampling time is used to determine the time step and the total number of 
measurements based on the simulation time. 
delta ? deltatime * 60 
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Number of Time Steps (nsamples) 


The number of time steps is determined using the value of delta and the total simulation 
time. This value is used by the main program loop for the total number of measurements. 
nsamples ? round {t! delta) 

Transition Matrix (F) 

The Transition Matrix is used to iterate target motion through time based on the time 

delta. The F-matrix is multiphed by the apphcable state matrix or truth state to determine 

the next target position [Ref. 9]. 

?1 delta 0 0 ? 

7 7 

AO 1 0 0 A 

7 I 

■ ?0 0 1 deltal 

7 7 

■X) 0 0 1 7 

Target Matrices (A) and (All 

The Target Matrices result from reading the TARDAT.MAT file and separating the 
targets into maneuvering (A) and non-maneuvering (Al) targets. This separation is 
accomphshed by performing a simple test: If the turn time in the last position of each 
target vector is a zero, the target is non-maneuvering regardless of the leg time input. 
The sensor platform and all targets (maneuvering and non-maneuvering) have the same 
format. The targets have emitters instead of sensors. For all non-maneuvering targets, 
the (F) matrix is used to transition the state vector. For all maneuvering targets, the 
(Fturn) transition matrix is used while the target is maneuvering. 
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ai1n/s, y ? pos, EIW,xl pos, heading, speed, altitude, type, name 


sensorl, sensor!, sensor!, sensor 4, comm#, legtime, L / Rturns, turntime 1 


Example: 


?o 

0 

0 

0 

270 

175 

25000 

1 

1 

1 

7 

9 

11 

0 

10 

1 

.43? 

•0 

A ? ? 

400000 

0 

300000 

170 

350 

25000 

1 

1 

1 

5 

6 

8 

10 

9 

1 

Asl 

?0 

o 

500000 

0 

200000 

240 

400 

20000 

1 

1 

1 

5 

6 

8 

10 

7 

0 

.5? 

o 

90 

200000 

0 

400000 

330 

100 

5000 

1 

1 

1 

1 

1 

1 

1 

8 

0 

.2 7 


?o 

250000 

0 

500000 

070 

500 

30000 

1 

1 

1 

1 

1 

1 

1 

8 

0 

0? 

A1 ? -0 

60000 

0 

60000 

355 

12 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0- 


700000 

0 

100000 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0? 


Sensor Platform Position (Sp) and (Spi) 

The sensor platform position is initialized with (Spi) by extracting the x-y position, the 
heading and the velocity from the first target in the TARDAT.MAT file. The sensor 
platform is always the first row vector of the TARDAT.MAT file. The sensor position 
used in all calculations of target motion is based on the navigation solution of the sensor 
platform (Sp). The sensor platform is sometimes referred to as Ownship. 


Ownship Sensors (ownsens) 

The sensors associated with the sensor platform (ownship) are determined based on user 
preference in the TARDAT.MAT file. Vector positions 10 - 13 contain the preferences 
for sensors as defined in the SENSERROR function. 

ownsens ? 1 7 9 11? (See Appendix A for a complete description of the sensors) 
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Maneuvering Target Transition Matrix (Fturn) and (ownfturn) 


The (Fturn) matrix is used to iterate target motion through time whenever the target or 
sensor platform is in a maneuvering condition. The matrix is set up by the FTURN 
function, which determines if the target is turning left or right and assigns the turn g’s. 
The (Fturn) matrix is based on each target’s heading, speed, and turn g values. For this 
simulation, the (Fturn) matrices are stacked for each target resulting in a 4*(# targets) x 4 
matrix. The (Fturn) matrix is used only for generation of truth data. 

||a|| 

? = — (turn rate in radians/sec); ? ? delta 

V 


?. sin(? ?) 

ol - 

• ? 

fturn = L'l {(si - A) }= e^' = '‘P I 

1 ? cos(? ?) 

? t 

k) sin(? ?) 


_ l?cos(??)? 

0 ^ ? 

0 ?sin(??)- 

^ sin(? ?) 7 

7 7 

0 cos(? ?) f 


where A is the System Matrix [Ref. 1,9]. 


Sensor Numbers (nsens) 

The (nsens) vector holds the values corresponding to the desired sensors of the sensor 
platform. The variable (ns) holds the total number of sensors associated with the sensor 
platform. The number of sensors can vary from 1 to 4. When more than one sensor is 
associated with the sensor platform, many of the vectors and matrices in the simulation 
are given a third dimension to allow analysis of the additional sensors. 
nsens ? 1 7 9? (See Appendix A for a complete description of the sensors) 
ns = 3 (in this case) 
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Tracker (trkr) 


The (trkr) vector is returned from the SETUPTRACKERS function containing the tracker 
that corresponds to each sensor. If the primary tracker is set to 1 by the user, the 
sensor 1/tracker 1 pair is set to 1. The 1 corresponds to the Probabihstic Data Association 
Eilter. Subsequent pairs are automatically set by the SETUPTRACKERS function, 
trkr = [sensorl/trackerl pair, sensor2/tracker2 pair, sensor3/tracker3 pair,...] 
trkr ? ?2 3 4 1?; 

In this case, sensor I corresponds to tracker 2 (IMM), sensor 2 corresponds to 
tracker 3 (Kalman), sensor 3 corresponds to tracker 4 (Const Gain Kalman), and sensor 4 
corresponds to tracker I (PDA). 

2. Initialization 

Initial x (xin) and (xim) 

The X vectors (xin) and (xim) are based on the initial values input by the user as the initial 
position, velocity and heading of the target. To keep them separate, xin contains the 
initial values for the non-maneuvering targets, and xim contains the maneuvering target 
data. The vectors contain the initial true states stacked for each target. The x2n and x2m 
vectors are simply the next truth data point generated by the TIMESTEP function. The 
TIMESTEP function simply multiphes the (E) matrix by the 4x1 truth state vector. 
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xin = 


1x1 
7 7 

'rfVX'rf 

lyl ’ 

7 7 

ivy) 


x im = 


1x1 
7 7 

7 7 
7VJ7 


x2n = F*xin = 


?x2? 

? .? 

? 3 ; 2 ? 
7 7 
7Vj27 


Maneuver Time Matrix (mantime) 

The Maneuver Time Matrix (mantkne) is set up by the MANMATRIX function and 
contains a column vector of times for each maneuvering target. The vector times are 
based on the leg times and turn times as defined by the user in the TARDAT.MAT file. 
For each target, the column vector holds the cumulative leg time and turn times. During 
matrix construction, if the max im um time is greater than the total simulation time, the last 
element is set to the maximum time. The size of the matrix is kept consistent by filling in 
the shorter vectors with zeros. 


7480 600 400 480 300 600? 

A580 750 500 540 450 700A 

mantime 1 • 

7900 900 900 900 750 900? 

7 0 0 0 0 900 0 7 


Initial State Estimation (xhatl) 

The initial state estimation for each target is determined by the INIT function and 
returned with the applicable sensor covariance (R) and prediction covariance (P). Also, 
the measurement error (derror) is returned for plotting and comparison purposes. The 
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additive plant noise term, Q is also returned. The INIT function is called for each sensor. 


??3 

?2 

2 

0 

7 

0 7 

7 

?1 

0 

0 

?? 

7 

0 

0 1 

?1 

0- 

0 

71 

? 2 

7 


?3 

7 

72 7’ 

• 9 

^■fo 

1 

7 

o' 

? 0 

0 

— 

—7 

7 

1 


? 


3 

2 7 

?0 


0 

jo 

0 

?2 

2 

7 ? 

■ ? 

7 

7 



The function POLE2CART is called within INIT for each target to obtain the 
noisy measurement (zcart), the sensor covariance (R), and the measurement distance 
error (disterror). Using the second set of truth data, the POLE2CART function is called 
again for use in Eease Squares Initiahzation. To simphfy initiahzation requirements, the 
least squares initialization is used for all trackers in the simulation, regardless of tracker 
type. The INIT function simply calls POEE2CART for each target using the appropriate 
truth position (xin/m), sensor error values (sigr, sigb, Sv) and sensor position (Sp). 
Within POEE2CART, the sensor position (sp) is subtracted from the target true x-y 
(ztrue). Next, the actual range and bearing (r,b) are calculated to the target and the 
measurement noise is added. The values are converted back to Cartesian coordinates and 
the sensor position is added back in to get the plotted position of the target. Einally, the 
measurement covariance in Cartesian coordinates is calculated and the distance error is 
calculated for determining the measurement error. 

z ? ztrue ? sp ; r ? Vz(l)" ? Zilf 

.?tani^| 

? z(V, ? 
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r = r + sigr*randn 
b = b + sigb*randn 

7rcos(b)7 lcos(b) ? rsm(b)l 

Z?7 9 sp; /v:?7 7 

9rsm(b)9 7sin(Z7) rcos(b) 7 

R = fx*Sv*fx’; ztilde = ztrue - z 

disterror ? V ztilde^.* ztilde 


Slight biases are introduced with this type of measurement processing [Ref. 1]. 
These errors are ignored in this simulation. Once 2 noisy measurements are obtained, the 


Kalman initiahzation can commence to create xhat and the covariance term: 


?v2? 


yflipped 


?z(2)? 

7 7 

7Z(1)7 


? .? 

■ ?vl? 
7 7 

7yl7 


? R2 zeroj'(2,2) ? 

COV ? 7 9, 9, 7 

izeros(2,2) Rll H * * Q* 

P = D*cov*D’ 
xhat =D*yflipped 


Sensor Platform Initial State (Spst) 

The initial state of the sensor platform is determined with the SPINIT function similar to 
the initiahzations of the targets. The sensor position state is determined using the same 
methodology as xhat for the targets. 
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3. Truth Data 


Non-Maneuvering Target Time Step 

If the target is not maneuvering at any time during the simulation, generation of the truth 
data is a simple step: 
xnew = F*xold 


Maneuvering Target Time Step 

For maneuvering targets, the truth data generation is a httle more comphcated. First, a 
test must be applied to see if each target is in a maneuver condition. This is where the 
previously described (mantime) matrix is used. If the target is in a maneuver state at the 
current simulation time: 


Target 5 


?480 

600 

400 

480 

300 

600? 


A580 

mantime ? • 

750 

500 

540 

450 

lool 

^ _ Turning 

Condition 

?900 

o 

900 

900 

900 

750 

900? 

o 

? 0 

0 

0 

0 

900 

0 7 




simt = 400 

xnew = ftum*xold 

When not in a turning condition, the previous expression is used for straight line motion. 


Sensor Platform Time Step 

The Sensor Platform is treated separately from the targets. To generate the truth data, the 
previous truth data is sent to the SPITMEAS function along with maneuver times and the 
ownship turning motion matrix (ownfturn). 
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spxtrue = F*spx2; or, spxtrue = ownfturn*spx2 (if maneuvering) 


4. Measurement 

Once all the truth data is generated, the measurement phase can commence. The 
truth data is sent to the POLE2CART function for each sensor and for each target to 
create the noisy measurements. The POLE2CART function calculates the measurements 
exactly the same as the measurements in the initialization section. 


Measurement Truth Data (zre) 

The measurement truth data is simply the ztrue vector reshaped into a stack of x-y row 
vectors. This reshaping is done to make the data more convenient to use with 
POEE2CART. 

7x1 yl7 


^x2 yl^ 

7x3 y37 
zre ? ? 7 

? • • ? 

? ? 
? ? 
Hxn yn7 


Sensor Platform Position Update (Sp) 

One of the key elements in the POEE2CART function is the sensor platform position 
because aU errors are measured with respect to the position of the sensor. To simulate the 
sensor platform navigation solution, the truth data is treated the same as the targets. Eirst, 
the truth data is sent to POEE2CART with some r, ? errors to create the measurement 


45 



noise (in the SPITMEAS function). This simulates errors in the platform’s GPS/Inertial 
navigation system. Next, the noisy measurement is sent to a Kalman Filter tracking 
function to create the navigation state update. 

Sensor Position State (spstate) 

The sensor position state is updated by calhng the KALMANl function with the previous 
state (Spst), the prediction covariance (spP), the measurement covariance (spR), the noise 
term (spQ), and the current noisy measurement (spmeas). The state prediction and 
update equations follow: 

Kalman Prediction: 

Pnext = F*spP*F’ + spQ 
xnext = F*xprevious 
Kalman Update: 

K = spP*H’*(H*spP*H’ + spR)'^ 

spP = (I - K*H)*Pnext*(FK*H)’ + K*spR*K’ 

spstate = xnext + K*(spmeas - H*xnext) 

Sp = H*spstate 

The (Sp) variable is the estimated x-y position as extracted from the state vector. 

Noisy Target Measurements (Zl) 

For each sensor (with associated errors), a call to the MEAS function results in a vector 
of noisy measurements for each target. The MEAS function simply calls POEE2CART 
for each target using the appropriate truth position (zre), sensor error values (sigr, sigb. 
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Sv) and sensor position (Sp). Within POLE2CART, the sensor position is subtracted 
from the target true x-y (ztrue). Next, the actual range and bearing (r,b) are calculated to 
the target and the measurement noise is added. The values are converted back to 
Cartesian coordinates and the sensor position is added back in to get the plotted position 
of the target. Finally, the measurement covariance in Cartesian coordinates is calculated 
and the distance error is calculated for determining the measurement error. 

For each target, the measurements, measurement covariance, and distance errors 
are stored (zret, R, derror) and returned to the for further storage and manipulation. 


?vl x2 x3 . 
7 yl y2 y3 . 



? ? 


5. Tracking 

The next step in the program is the tracking section. Four different trackers are 
included in the current configuration of the fusim program. A complete hsting of the 
MATFAB code for all the algorithms is contained in Appendix B. The TRACKER 
function is called for each sensor with the (trkr) vector for determining which tracker is 
associated with which sensor. The TRACKER function makes the call to the appropriate 
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tracking function. The tracking algorithms used mfusim are hsted in order of precedence 
starting with the Probabahstic Data Association Filter. 

Probabilistic Data Association Filter (PDAF) 

The PDAF [b-s] is the only tracking code included with fusim that allows analysis of 
targets moving in clutter. This implementation of the PDAF is not designed for use with 
multiple targets or maneuvering targets. However, increasing the additive noise term 
(Qp) allows sufficient tracking through turns for this analysis. 

Initially, the probabihty of detection (Pd) and the gate density probabihty (Pg) are 
set. The gate probability is usually around 1 [Ref. 1]. The number of targets and clutter 
points (m) are also predetermined. For this PDAF, one target is always assumed to be 
present. Another simphfication for this filter is the calculation of the tracking gate. For 
each clutter point within the gate, a simple call to POLE2CART is used to generate the 
point. This assures that the clutter point always falls within the gate. These values do not 
change in this adaptation of the Probabilistic Data Association Filter. 

Pd = 1.0; Pg = .99; m = 2 

For the gate calculation: 

? range ? 900 fcCt 

• bearing ? 2 deg 

, ?? \ange 0 ? 

gate variance '9 ^ ? 

7 0 ? hearing 7 

qsquared = 1000 
Qp = qsquared*Q 
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The weighted predictions: 

X Fx 

./V J • -/ 

Pnext ? FPpF^il Qp 

For calculation of the gate volume factor, the semi-major and semi-minor axis in feet: 
minor axis, a = 600; 
major axis, b = 1500; 

Vk ? 2?ab 


7 7 — 

■ Vk 

PreS ? HPnextP'. 

Calculate the innovation covariance (Sk) for the actual target, based on the measurement: 
Sk ? PxeS ? R 
K ? PnextH'Sk^^ Gain 


bl 


7, |iV7l?PdF>a7 

?? 2?5k 2 p— —ICA7 

7 ' ' hi Pd 7 


sum = b 
? ? zret ? Hx^ 

?f? ^Sk■'^ ? 

5Mm ? 5Mm ? el 

For each clutter point, the POLE2CART function is called to create the noisy 
measurements and the association probabihties (?) are calculated. 
h j h zret ? //v, 
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sum ? sum ? Cj 

6: 

7 , 7 _ I— 

^ sum 

Then, for the combined innovation calculation: 

(y ry ^ ry ry 

' k ' • / ? 1 ’ * ’ * 

For determining the spread of the innovations term, the calculations are broken down as 
follows: 

For each clutter point (and target), 
spread ■>‘7 ?,?,????,?? 

P ? K^lspread’lKl 
Pc ? Pnext ? KSkK'^ 

7 7 ^ 

• 0 • 

sum 

Pp ? iJiPclP 

Final state estimation: 

xl x,lKl, 


Interacting Multiple Models 

The 2-model Interacting Multiple Models filter allows tracking of a maneuvering target 

based on a straight hue motion model and a turning model. Initially, the probabihty of 
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conducting a turn if in straight line motion (alpha) and the probability of stopping a turn 
(beta) are set. This allows setting up the Markov 2x2? matrix. 

71 ? ? 7 7 

? ? 7 7 

7 ? 1? ?7 

In this tracker, ? =0.1 and ? = 0.33333. These values can be changed for user 
preference. 

Initial State Likelihood: (straight hue track) 

? 1 ? 

? ? 7 7 ; Initially, ? i = 1.0, indicating the target is assumed to be in straight line motion 
tO? 

(SLM). The state li kelihoods are recalculated each time the tracker is called for a new 
update. 

Set the noise terms for the two different models: 

Qi = L0*Q; 

Qturn = 10000*Q; 

Pre-process the cbars and mu's: 


7 7 ^ cy cy ^ 

Cj . . JJ . J . .21-2 


The initial xhat guesses are based on the probabihties that the target will continue 
straight hue if in SLM, turn if SLM, continue turning if in turning motion, or return to 
SLM if turning: 

7777 7777 

? Cl ? ? Cj 7 
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7 7 7 

7 7 7 7 7 


9 C 7 • 12 • 1 


^ z 2 

• ^17 _ 

7 • ^2 7 _ 7 


? ^2 

? ? C2 7 


7 7 

7 7 

T 

9 • 11 • 1 . 

9- 9 • 21 • 2 

* 11 

• 9 

• 21 • — 


Cl 

^1 


7 7 

? ? 

T 

9 • 12 • 1 . 

9- 9 • 22 • 2 

• 12 

• 9 

• 22 • — 


C2 

^2 

Vii 

X "> X • 

• -^1 • 5 

V21 ? V2 ? 

^12 

X "> X • 

• -^1 • -^Z2 ’ 

V22 ? X2 ? V^2 


p ^ X ^ Pti/f m 77~ ^P ^ X ^ Ptu m ^ x jv ^ ^ 

-*01 • • 11 •-* • '^ ll'^ll • • • 21 •-* • '^ 21'^21 •’ -* 02 • • 12 •-* • '^ 12'^12 ’ ’ ’ 22 -* • '^ 22*^22 ’ 

Next, the weighted predictions are performed: 

/V r> T~<^ ^ O T~'^ 

Xi ? Fx^i; ^2 ? Fx^2 


Pi ? pPoi^"^? 

P2 ? FP^F^l Qturn 

Before proceeding to the next step, the straight hne Kalman gain, turning gain, and 
covariances are determined separately: 

K ? P^HtHP,HV. P?‘ 

Kturn ? P^HtHP^HV. rT 


c 

Pi 7 In KH ?Pi ?/ 7 KH? 7 KRK7 

c 

Ptum ? ?/ ? KturnH‘I p^’Ii 7 KturnHl' 7 KturnRKturn7 
5; ? HP,H\7 R 
7 HP^H\7 R 
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The measurement update step is conducted by calculating the state estimate for straight 


hne and turning, and then combining the results for plotting: 
Straight hne 
Zj ? zret ? Hx^ 

X X K? 

■^slm • -'•1 • ^^1 

Turning 

Z 2 ? zret ? Hx 2 

^,urn ? ^2 ? Kturnz 2 


Combining the result for plotting and error analysis, the state estimate: 

r 7 ? r 7 7 V 

• • 2-^tum 


Score the results (simplified with only two models): 


7 

• 1 


9 1 T 2 ’ 


2? 51 2 


2 ?r7c9?lr ?2 

Z 2 77 

2 ■ 2 9 22 

77 2 99 

e 


2 ? ||52||i 


C ‘t c c 

c . . [Cj . . 2^2 


7 C 

7 7 • 1^1 • 

• 1 • ’ 

C 


7 C 

7 7 •2^2 . 

• 2 • ’ 

C 


77 7 

? ?? ‘? 
?? 2 ? 


For the next prediction/update: 


o ^ 

-^1 ? ; 


^ o ^ 
-^2 ? 
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Kalman Filter 


The Kalman filter used for tracking the targets is exactly the same code as the equations 
used in the sensor platform navigation solution. 

Kalman Prediction: 

Pnext = F*spP*F’ + spQ 
xnext = F*xprevious 
Kalman Update: 

K = spP*H’*(H*spP*H’ + spR)'^ 

spP = (I- K*H)*Pnext*(FK*H)’ + K*spR*K’ 

spstate = xnext + K*(spmeas - H*xnext) 

Sp = H*spstate 

Constant Gain Kalman Filter 

The constant gain Kalman filter is the simplest case included in the. fusim program. 
qsquared = 1000.0; 

Q = qsquared*Q; 

First, the constant Kalman gain is determined using the MATLAB function DLQE for 
each target being tracked. This value of K is used for aU subsequent state calculations. 

X Fx 

■^pred ■ ^ 

^ ? ^pred ? Kizret ? 
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6. Data Link Simulation 


For the data link simulation, bias error and noise errors are added to the state of 
the primary sensor on the sensor platform. The errors are added to simulate reference 
frame differences, noise in another platform, and transmission errors. The link report is 
then formatted to allow many different types of data to be sent as part of the report, 
including attribute data, and the time of the report. For each target: 

X = X + randn*50 + 100 feet 
y = y + randn*50 + 100 feet 
vx = vx + randn*10 feet/sec 
vy = vy + randn*10 feet/sec 
linkvector = [x,vx,y,vy,hdg,time, . ] 

The hnk reports are meant to be used in a track-to-track data association 
algorithm, which is not available in this first edition of fusim. The x and y data are the 
only data used in this simulation for plotting purposes. 

7. Storage and Plotting 

For every instance of target motion, the truth data (posout), “measured” position 
(zout), tracker state (state), and least square error (Lserror) are stored in a matrix. The 
same data is collected for the sensor platform. After each Monte Carlo run, the current 
data are added to the previous cumulative data for eventual averaging over the number of 
runs. Three different plotting functions are provided for plotting the real-time results or 
for plotting the tracks and the errors for each target and sensor. 
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D. EVALUATION PARAMETERS 


1. The Baseline Target Set 

The target set used for all simulations are presented in Table 2. This setup was 
designed to be flexible enough to add or remove maneuvering or non-maneuvering 
targets via the TARDAT.MAT file. The format of the TARDAT.MAT file must be 
followed precisely, with zeros entered where no data is desired. 


Track 

Initial Position 

(x,y) 

Initial 

Course 

(degrees) 

Initial 

Speed 

(knots) 

Maneuvering? 

(Y/N) 

Notes 

Ownship 

0000/0000 

270 

175 

Y 

Airborne surveillance platform 

Target 1 

400000/300000 

170 

350 

Y 

Maneuvering Fighter Aircraft 

Target 2 

500000/200000 

240 

400 

Y 

Maneuvering Fighter Aircraft 

Target 3 

200000/400000 

330 

100 

Y 

Slow-Moving Air Track 

Target 4 

250000/500000 

070 

500 

N 

Air Liner, Non-Maneuvering 

Target 5 

60000/60000 

355 

12 

N 

Ship 

Target 6 

700000/100000 

0 

0 

N 

Fixed Target/Emitter 


Table 2. 

The target set used 1 

for all simulations. 


2. Simulation Set 

For this initial edition of fusim, several different sets of runs were performed to 
examine the vahdity of the code, to check the simulation in a multiple sensor/multiple 
target scenario, and to observe algorithm behaviors in a simple combined solution. The 
three major tests are described as follows: 

Test 1: Track algorithm vahdity - The vahdity of the simulation was checked by 
running single sensor/multiple track test sets and comparing the mean measurement noise 
to the tracking solution for each track. Initially, zero error was input to check for 
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problems in the code. For the next set of runs, the noise term (q ) was adjusted in two of 
the trackers to compare straight line motion (SLM) and maneuvering target tracking 
characteristics. Finally, the track solutions were compared by running three different sets 
of runs using a single sensor with all four trackers. Table 3 shows a breakdown of the 
tests and the test subsets. 

Test 2: Multiple Sensor/Multiple Target Tracking - Test 2 was run to examine 
the use of multiple sensors and tracker types in a multiple target tracking scenario. 
Representative plots were produced to show the differences in mean measurement noise 
and tracking solutions when using four different sensors and four different tracking 
algorithms. This test was simply a demonstration of the multiple sensor/multiple target 
case and a precursor to the next test. 

Test 3: Simple Combined Solution - A third group of test sets were run to 
demonstrate a simple combined tracking solution using the multiple sensor/multiple 
target tracking results of Test 2. The combined state result for each track was simply a 
weighted average of each of the individual tracker states. The weights were selected 
empirically based on the tracker performance results observed in Test 1. 

The combined state estimate was based on the following weights (percentages) 
for each of the sensors: 


si: 

Sensor 1 = 0.3 

(Surveillance Radar) 

s2: 

Sensor 2 = 0.4 

(IFF System) 

s3: 

Sensor3 = 0.1 

(ESM System) 

s4: 

Sensor 4 = 0.2 

(IRST System) 


Then, for the combined state estimate: 


xhat = sl*xhat(sensorl) + s2*xhat(sensor2) + s3*xhat(sensor3) + s4*xhat(sensor4). 
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For all of the tests, the runs were conducted with the same target set and target 


initial parameters (Table 2). Table 3 is a description of the evaluation sets. The sensor 
descriptions and associated errors are included in Table 4. 


TEST 

NAME 

TEST SUBSET 

Test 1 

Track Algorithm 
Validity 

Set (a) Zero sensor error 

Set (b) Analyze algorithm noise terms 

Set (c) Airborne surveillance radar 

Set (d) Surface fire control radar 

Set (e) Infrared search and track system simulator 

Test 2 

Multiple Sensor 
Multiple Track 

Set (a) Four different sensors and trackers 

Test 3 

Observations on 
Algorithm 
Behavior in a 
Simple Combined 
Solution 

Set (a) Single sensor type/combined solution with four 
tracking solutions 

Set (b) Four sensors/combined solution 

Set (c) Four sensors/combined solution, modified noise 
terms for turn performance 

Tal 

ble 3. Fusim Evaluation Set. 
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Sensor 

Range 

Error 

Bearing 

Error 

Notes 

Airborne Surveillance 
Radar 

100 ft 

.005 rad 

Simphfied model of an airborne 2D 
radar. 

IFF Detection System 

50 ft 

.001 rad 

Simple interrogation and reply 
range and bearing only. 

ESM System 

500 ft 

.01 rad 

A simphfied model of an ESM 
system. 

IR Search and Track 

6000 ft 

.0001 rad 

A simphfied model of an IRST 
system. 

Surface Fire Control 
Radar 

10 ft 

.0001 rad 

A simplified model of a Fire 

Control Radar system. 


Table 4. Test Sensor Descriptions. 


3. Test Set-up Parameters 


(a) 15 minute simulation time. 

(b) 100 runs (Monte Carlo Simulations). 


Unless otherwise noted, the simulation parameters of Tables 2-4 are used throughout the 
simulation. 
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V. SIMULATION RESULTS 


The plan position indicator display of the tracking results of Table 2 is included in 
Figure 2. This is a representative plot showing the location of the targets and the sensor 
platform. The same figure was produced for every test case that was run, but was not 
included as a figure for every case. Unless otherwise noted, the targets referred to in the 
rest of this section are the targets indicated in Figure 2. 


X 10* Mean State (Tracking Solution) for All Sensors(20 runs) 



LONGITUDE - in feet for now ^ .jg* 


Figure 2. Plan Position Indicator (PPI) of Simulation Results. 
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A. TEST 1: TRACK ALGORITHM VALIDITY 

Track Algorithm Validity Set (a) Results 

To prove that the simulation was working for generation of truth data, 
measurement data and tracking data, the code was initially run with zero error input for 
the range and bearing errors. The sensor platform had zero error, the sensors each had 
zero error, and the tracker was allowed to produce a track and determine the Least Square 
Error compared to the truth data. With zero error in the sensor, several runs were 
conducted to demonstrate stabihty in the tracking codes. Under ideal conditions, with no 
error in the sensors, the track produced had zero error for straight line motion (Figure 3). 
Errors greater than zero were noted for maneuvering targets, but only in the turn. As 
expected, the errors quickly returned to zero upon completion of the turn. 

Track Algorithm Vahdity Set (b) Results 

For the next set of runs, the errors were reinstated and the noise terms in each of 
the trackers were adjusted to observe the differences in the tracking solution. For the 
PDAF, one clutter point was used in the algorithm. The number of clutter points was 
fixed at one to simplify the code and force the tracker to deal with a clutter point at every 
iteration. This would not be the case in an actual PDAF, where normally a test is 
performed at every measurement update to see if a clutter return originates from within 
the tracking gate. For the IMM, the straight line motion (SFM) model and turn model 
noise terms were adjusted to find a solution that performed well in the turn. The Kalman 
Filter (KF) and the Constant Gain Kalman Filter (CGKF) were not evaluated during this 
portion of the test. 
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Test 1: Set (a) 

Tracking Algorithm Validity 

Description: 

Determine if the tracking algorithms are valid with zero tracking error 

Sensor 

Airborne Surveillance 
Radar 

Identification Friend or 

Foe System 

Electronic Surveillance 
System 

Infrared Search 
and Track 
System 

Tracker 

PDAF 

IMM 

Kalman Filter 

CGKF 






Settings: 

1000.0 

-SLMq^= 1.0 

- q"' = 1000.0 

-q^= 1000.0 


- 1 clutter pt 

- Turn = 10000.0 



Targets 

Maneuvering Aircraft 

Non-Maneuvering Aircraft 

Slow-Moving Ship 

Fixed Site 


Note: The test conditions are highlighted in each figure test table. 



IQ''* Mean Error (-) vs. Time, 1 Sensor Type/4 Trackers/1 run 

1 

- Zero added error 
-All Trackers 


Non-maneuvering 

Target 

Error, (feet) 

o 


-1 

I _1_1_1_1_1_1_1_1_1_1 

354.14354.15354.16354.17354.18354.19 354.2 354.21354.22354.23354.24 

Time, (seconds) 


Figure 3. Tracking Solutions with no Measurement Errors. 
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The airborne surveillance radar was used to check the vahdity of the PDAF. The 
PDAF was isolated as a tracker, the noise term was adjusted to intentionally force the 
tracker to miss the turn, and the results were plotted to demonstrate the subtle weakness 
of the PDAF. The PDAF worked well for tracks in straight-hne motion with q =1.0, but 
lost the track in the turn (Figures 4 and 5). Instead of updating the track state with the 
measurement, the track state was updated with a more heavily weighted clutter point. At 
q = 100 (Figure 6), the PDAF was able to maintain a reasonable track through the turn 
for this set of runs. However, for higher accelerations or turn rates, this PDAF setup 
could not be expected to maintain the track. In fact, the model runs a rather low turn 
acceleration (1.5g) that does not represent a worst-case scenario. For this simple 
implementation of the PDAF, it is always possible for the track to get pulled off by the 
clutter point, mainly because of the previous assumption: The clutter point always falls 
within the tracking gate (recall that the tracking gate is centered around the predicted 
position of the target). 

The mean measurement noise appears to be dropping off rapidly in each of the 
figures. This is due to the change in geometry from the moving target and sensor 
platform. At q = 250 (Figure 7), the PDAF tracked fine through the turn, but the mean 
measurement error was about 400 feet higher for the same target. These adjustments for 
the noise term are the same as the noise adjustments in the Kalman filter. For the next set 
of runs, the noise term was adjusted to allow the PDAF to track through turns in search of 
an optimum solution. At q = 250, the average track error was up around 700 feet for 
SLM, but down around 300 feet with q =1. To ensure a constant turn performance 
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throughout the rest of the tests, the PDAF was run with a noise term of q = 1000 (unless 
otherwise noted). 

For the IMM, the two models were varied by simply altering the noise term in 
each model. For the straight hne motion model, the noise term was very small to take 
advantage of the SLM performance of the Kalman filter. The Markov probabilities were 
0.1 for the turn probabihty, and 0.33 for the stop turn probability. Changing these 
probabihties do affect the performance of the filter, but this aspect of the IMM was not 
rigorously examined in this thesis. The IMM was run for several different cases with 
varying values of the noise term. For the turn model, the noise term was gradually 
increased until the filter performed consistently well against a maneuvering target. Only 
two of the cases are displayed in Figures 8 and 9. Reasonable performance was obtained 
in the turn with the SLM and turn q = 1 and 10000 respectively. 

The noise terms in the Constant Gain Kalman Filter (CGKF) and the Kalman 
Filter (KF) were not evaluated for the purpose of this Thesis due to proven performance 
in other exercises [Ref. 9]. 


65 



Testl: Set{b) 

Tracking Algorithm Validity 

Description: 

Determine an applicable noise term (q^) for follow-on tests 

Demonstrate susceptibility to losing the track in a turn of the PDAF with a low noise term 

Sensor 

Airborne Surveillance 
Radar 

Identification Friend or 

Foe System 

Electronic Surveillance 
System 

Infrared Search 
and Track 
System 

Tracker 

Settings: 

PDAF 

IMM 

Kalman Filter 

CGKF 

1.0 

- 1 clutter pt 

-SLMq^= 1.0 

-Turnq^= 10000.0 

- q- = 1000.0 

- q- = 1000.0 

Targets 

Maneuvering Aircraft (2) 

Non-Maneuvering Aircraft 

Slow-Moving Ship 

Fixed Site 


"IQ® Mean State (Tracking Solution) for 1 Sensor/PDAF(20 runs) 



Figure 4. Tracking Results for PDAF with Low Noise Term. 
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Testl: Set{b) 

Tracking Algorithm Validity 

Description: 

Determine an applicable noise term (q^) for follow-on tests 

Demonstrate susceptibility to losing the track in a turn of the PDAF with a low noise term 

Sensor 

Airborne Surveillance 
Radar 

Identification Friend or 

Foe System 

Electronic Surveillance 
System 

Infrared Search 
and Track 
System 

Tracker 

Settings: 

PDAF 

IMM 

Kalman Filter 

CGKF 

-q^= 1.0 

- 1 clutter pt 

-SLMq^= 1.0 

-Turnq^= 10000.0 

- q- = 1000.0 

- q- = 1000.0 

Targets 

Maneuvering Aircraft 

Non-Maneuvering Aircraft 

Slow-Moving Ship 

Fixed Site 


6000 r 


5000 - 


4000 - 


O) 

Q> 


3000 - 


Mean Error (-) vs. Time, 1 Sensor/PDAF/20 runs 


Mean Measurement Error 


PDA tracked off 
with the clutter point 


2000 


1000 I 


\ 

V/i'i 


I_I_I_I_I 


0 100 200 


300 400 500 600 

Time, (seconds) 


700 800 900 


Figures. PD AF with Low Noise Term. 
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Testl: Set{b) 

Tracking Algorithm Validity 

Description: 

Determine an applicable noise term (q^) for follow-on tests 

Demonstrate how increasing the noise term improves the turn performance of the PDAF 

Sensor 

Airborne Surveillance 
Radar 

Identification Friend or 

Foe System 

Electronic Surveillance 
System 

Infrared Search 
and Track 
System 

Tracker 

Settings: 

PDAF 

IMM 

Kalman Filter 

CGKF 

100.0 

- 1 clutter pt 

-SLMq^= 1.0 

-Turnq^= 10000.0 

- q- = 1000.0 

- q- = 1000.0 

Targets 

Maneuvering Aircraft 

Non-Maneuvering Aircraft 

Slow-Moving Ship 

Fixed Site 



Figure 6. PDAF with Increased Noise Term. 
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Testl: Set(b) 

Tracking Algorithm Validity 

Description: 

Determine an applicable noise term (q^) for follow-on tests 

Demonstrate how increasing the noise term improves the turn performance of the PDAF 

Sensor 

Airborne Surveillance 
Radar 

Identification Friend or 

Foe System 

Electronic Surveillance 
System 

Infrared Search 
and Track 
System 

Tracker 

Settings: 

PDAF 

IMM 

Kalman Filter 

CGKF 

- = 250.0 

- 1 clutter pt 

-SLMq^= 1.0 

-Turnq^= 10000.0 

- q- = 1000.0 

- q- = 1000.0 

Targets 

Maneuvering Aircraft 

Non-Maneuvering Aircraft 

Slow-Moving Ship 

Fixed Site 


Mean Error (-) vs. Time. 1 Sensor/PDAF/20 runs 



Figure 7. PDAF with Increased Noise Term. 
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Testl: Set(b) 

Tracking Algorithm Validity 



Description: 

Determine an applicable noise term (q^) for follow-on tests 

Demonstrate how changes in the models affects the IMM performance. 


Sensor 

Airborne Surveillance 

Identification Friend or 

Electronic Surveillance 

Infrared Search 


Radar 

Foe System 

System 

and Track 
System 

Tracker 

PDAF 

IMM 

Kalman Filter 

CGKF 

Settings: 

- q" = 250.0 

-1 clutter pt 

-SLMq^= 1.0 
-Turnq’= 10000.0 

- q- = 1000.0 

- q- = 1000.0 

Targets 

Maneuvering Aircraft 

Non-Maneuvering Aircraft 

Slow-Moving Ship 

Fixed Site 



Figure 8. IMM Tracking Performance with Varying Models. 
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Testl: Set{b) 

Tracking Algorithm Validity 

Description: 

Determine an applicable noise term (q^) for follow-on tests 

Demonstrate how changes in the models affects the IMM performance. 

Sensor 

Airborne Surveillance 
Radar 

Identification Friend or 

Foe System 

Electronic Surveillance 
System 

Infrared Search 
and Track 
System 

Tracker 

Settings: 

PDAF 

IMM 

Kalman Filter 

CGKF 

- = 1000.0 

-1 clutter pt 

-SLMq^= 1.0 
-Turnq’= 10000.0 

- q- = 1000.0 

- q- = 1000.0 

Targets 

Maneuvering Aircraft 

Non-Maneuvering Aircraft 

Slow-Moving Ship 

Fixed Site 



Figure 9. IMM Tracking Performance with Varying Models. 
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Track Algorithm Validity Set (c) Results 


The airborne surveillance radar was the first sensor to be evaluated with all 4 
trackers (Figures 10 through 13). For the maneuvering target case of Figure 10, the 
trackers all performed reasonably well. However, the CGKF took far too long to 
estabhsh the track The large initial error (>6000ft in this case) is a characteristic of the 
CGKF. The only tracker with an advantage in this case was the IMM. It was shghtly 
worse in SLM and shghtly better in the turn. It would be easy to force the KF as the 
optimum solution in SLM by reducing the noise term to 1, but turn performance would 
suffer greatly. Also, the IMM could be forced to perform better in the turn by increasing 
the turn noise term. However, the SLM performance would be further degraded. 

In the case of the non-maneuvering target (Figure II), the PDAF and the Kalman 
filter were about the same. The IMM was shghtly worse mainly because of the 
expectation of a turn. The CGKF performed adequately, but errors in the CGKF 
increased much faster than the other trackers as range to the target and mean 
measurement error increased. The mean measurement errors increased over time due to 
increasing range from the sensor platform. 

In the case of the slow-moving ship (12 knots), the PDAF immediately picked up 
on the clutter point and tracked off (Figure 12). Some parameters may need to be 
changed in the PDAF for tracking slow moving targets. The PDAF may not be weh 
suited for tracking slow moving and stationary targets in a clutter region without 
modifications to make the code more adaptable. 
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For the fixed site (Figure 13), the trackers aU seemed to perform about the same 
for the surveillance radar, even though the mean measurement error was considerably 
higher. The IMM error was shghtly higher throughout, but that may stiU be due to the 
turn expectation and the high q of the turn model. 
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Test 1: Set (c) 

Tracking Algorithm Validity 

Description: 

Observe the performance of the trackers using the same sensor type with each of the trackers 

Sensor 

Airborne Surveillance 
Radar 

Identification Friend or 

Foe System 

Electronic Surveillance 
System 

Infrared Search 
and Track 
System 

Tracker 

PDAF 

IMM 

Kalman Filter 

CGKF 






Settings: 

1000.0 

-SLMq^= 1.0 

- q"' = 1000.0 

-q^= 1000.0 


- 1 clutter pt 

-Turnq’= 10000.0 



Targets 

Maneuvering Aircraft 

Non-Maneuvering Aircraft 

Slow-Moving Ship 

Fixed Site 


Mean Error (-) vs. Time, 1 Sensor Type/4 Trackers/lOO run(s) 



6000 r 


5000 - 


AIRBORNE SURVEILLANCE RADAR 
SIMUATOR 


1000 


0 100 200 300 400 500 600 700 800 900 

Time, (seconds) 


Maneuvering Target 


4000 


T3000 


2000 


Mean Measurement 
Error 


Figure 10. Tracker Comparison with One Sensor Type for the Maneuvering Target. 
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Test 1: Set (c) 

Tracking Algorithm Validity 

Description: 

Observe the performance of the trackers using the same sensor type with each of the trackers 

Sensor 

Airborne Surveillance 
Radar 

Identification Friend or 

Foe System 

Electronic Surveillance 
System 

Infrared Search 
and Track 
System 

Tracker 

PDAF 

IMM 

Kalman Filter 

CGKF 






Settings: 

1000.0 

-SLMq^= 1.0 

-q^= 1000.0 

-q^= 1000.0 


- 1 clutter pt 

- Turn = 10000.0 



Targets 

Maneuvering Aircraft 

Non-Maneuvering Aircraft 

Slow-Moving Ship 

Fixed Site 


6000 r 


5000 


4000 


Q> 


t 3000 


Mean Error (-) vs. Time, 1 Sensor Type/4 Trackers/lOO run(s) 

AIRBORNE SURVEILLANCE 
RADAR 


Mean Measurement 
Error 


\ 


Non-Maneuvering 

Target 

CGKF 



Figure 11. Tracker Comparison with One Sensor Type for the Non-maneuvering Target. 
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Test 1: Set (c) 

Tracking Algorithm Validity 

Description: 

Observe the performance of the trackers using the same sensor type with each of the trackers 

Sensor 

Airborne Surveillance 
Radar 

Identification Friend or 

Foe System 

Electronic Surveillance 
System 

Infrared Search 
and Track 
System 

Tracker 

PDAF 

IMM 

Kalman Filter 

CGKF 






Settings: 

1000.0 

-SLMq^= 1.0 

- q"' = 1000.0 

-q^= 1000.0 


- 1 clutter pt 

-Turnq’= 10000.0 



Targets 

Maneuvering Aircraft 

Non-Maneuvering Aircraft 

Slow-Moving Ship 

Fixed Site 



Figure 12. Tracker Comparison with One Sensor Type for the Slow-moving Ship. 
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Test 1: Set (c) 

Tracking Algorithm Validity 

Description: 

Observe the performance of the trackers using the same sensor type with each of the trackers 

Sensor 

Airborne Surveillance 
Radar 

Identification Friend or 

Foe System 

Electronic Surveillance 
System 

Infrared Search 
and Track 
System 

Tracker 

PDAF 

IMM 

Kalman Filter 

CGKF 






Settings: 

1000.0 

-SLMq^= 1.0 

- q"* = 1000.0 

-q^= 1000.0 


- 1 clutter pt 

-Turnq’= 10000.0 



Targets 

Maneuvering Aircraft 

Non-Maneuvering Aircraft 

Slow-Moving Ship 

Fixed Site 


Mean Error (-) vs. Time, 1 Sensor Type/4 Tracl<ers/200 run(s) 

6000 


5000 


4000 


4 3000 

o 

m 

2000 

1000 


0 

0 100 200 300 400 500 600 700 800 900 

Time, (seconds) 


Figure 13. Tracker Comparison With One Sensor Type for the Fixed Site Target. 


AIRBORNE SURVEILLANCE 
RADAR 


Mean Measurement 
Error 
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Track Algorithm Validity Set (d) Results 


For Set (d), the surface fire control radar was examined to look for anomalies in 
the trackers when paired with a highly accurate sensor (Figure 14). The ranges and 
accuracy’s of the sensor may be unreahstic, but there was no apparent advantage to any 
of the trackers with a sensor this accurate. In other words, the simple trackers worked 
just as well as the more complex trackers. This phenomenon gives the appearance that 
the state prediction and estimation was a waste of time and computing power. However, 
estimation is stiU necessary for predicting the state of the target, for tracking a 
maneuvering target, and for guiding a missile or projectile to an impact point. For the 
purposes of this test, the low errors were used to look for anomahes in the trackers. 
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Testl: Set(d) 

Tracking Algorithm Validity 

Description: 

Observe the performance of the trackers using the same sensor type with each of the trackers 

Sensor 

Fire Control Radar 
Simulator 

Identification Friend or 

Foe System 

Electronic Surveillance 
System 

Infrared Search 
and Track 
System 

Tracker 

PDAF 

IMM 

Kalman Filter 

CGKF 






Settings: 

1000.0 

-SLMq^= 1.0 

- q^= 1000.0 

-q^= 1000.0 


- 1 clutter pt 

-Turnq^ = 10000.0 



Targets 

Maneuvering Aircraft 

Non-Maneuvering Ahcraft 

Slow-Moving Ship 

Fixed Site 


Mean Error (-) vs. Time, 1 Sensor Type/4 Trackers/100 run(s) 

20 


18 


16 


14 


_ 12 

* 0 ) 

CD 

10 

o 

m g 


6 


4 


2 


0 

0 100 200 300 400 500 600 700 800 900 

Time, (seconds) 


FIRE CONTROL RADAR 


Mean Measurement Noise 



Trackers 


J_I_L 


Figure 14. Tracker Comparison with One Sensor Type for the Maneuvering Aircraft. 
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Track Algorithm Validity Set (e) Results 


The last test case in the first set of runs was the Infrared Search and Track (IRST) 
simulator. In formulating this very simple model of the an IRST system, the bearing to 
the target was assumed to be very accurate and the range was assumed to be very 
ambiguous. In a true IRST system, the range is not measured by the system, unless by 
some other means. To keep this analysis simple, it is assumed that the IRST has already 
estabhshed a track, and the measurement is based on the position report from the tracker. 

For the maneuvering target (Figure 15), the higher mean measurement noise of 
the IRST system made httle difference in the tracking performance. The CGKF error 
steadily increased as range to the target increased. An error increase in the other trackers 
was imperceptible. The IMM handled the turn fairly well, yet the Kalman and the PDAF 
both did worse in the turn. 
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Test 1: Set (e) 

Tracking Algorithm Validity 

Description: 

Observe the performance of the trackers using the same sensor type with each of the trackers 

Sensor 

Airborne Surveillance 
Radar 

Identification Friend or 

Foe System 

Electronic Surveillance 
System 

Infrared Search 
and Track 
System 

Tracker 

PDAF 

IMM 

Kalman Filter 

CGKF 






Settings: 

1000.0 

-SLMq^= 1.0 

- q"' = 1000.0 

-q^= 1000.0 


- 1 clutter pt 

-Turnq-= 10000.0 



Targets 

Maneuvering Aircraft 

Non-Maneuvering Aircraft 

Slow-Moving Ship 

Fixed Site 


8000 r 

7000 - 

6000 - 

5000 


Mean Error (-) vs. Time, 1 Sensor Type/4 Trackers/100 run(s) 

INFRARED SEARCH AND 
TRACK SIMULATOR 


Mean Measurement 
Error . 

\ 


Maneuvering Target 



100 200 300 400 500 600 700 800 900 

Time, (seconds) 


Figure 15. Tracker Comparison with One Sensor Type for the Maneuvering Aircraft. 
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For the non-maneuvering target (Figure 16), the results were similar to the 
maneuvering target case. The CGKF error gradually increased, while the other trackers 
remained steady around 1000 feet of error. For the slow moving ship of Figure 17, the 
CGKF error seemed to vary with distance from the sensor platform, even though the 
measurement error remained about the same over the simulation period. The PDAF, the 
IMM, and the KF aU hovered around 1000 feet of error. 

In this case of the fixed site (Figure 18), the PDAF lost the track, yet the other 
trackers did consistently well. The clutter point may have caused sufficient variation in 
the state estimate to induce a velocity and cause the solution to track off. This problem 
did not occur every time with the PDAF and may indicate a higher degree of uncertainty 
in a PDAF solution. To handle a slow-moving or non-moving target in a clutter region, 
some modification within the PDAF would be required to help classify the target and 
allow better estimates. For example, a velocity threshold could be used to help decide the 
class or type of target, so that unrealistic velocity estimates could be e lim inated from the 
solution. 
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Test 1: Set (e) 

Tracking Algorithm Validity 

Description: 

Observe the performance of the trackers using the same sensor type with each of the trackers 

Sensor 

Airborne Surveillance 
Radar 

Identification Friend or 

Foe System 

Electronic Surveillance 
System 

Infrared Search 
and Track 
System 

Tracker 

PDAF 

IMM 

Kalman Filter 

CGKF 






Settings: 

1000.0 

-SLMq^= 1.0 

- q"' = 1000.0 

-q^= 1000.0 


- 1 clutter pt 

- Turn = 10000.0 



Targets 

Maneuvering Aircraft 

Non-Maneuvering Aircraft 

Slow-Moving Ship 

Fixed Site 


Mean Error (-) vs. Time, 1 Sensor Type/4 Trackers/lOO run(s) 

8000 


7000 


6000 


5000 

"53 

Q> 

^4000 

o 

lii 

3000 


2000 


1000 


0 

0 100 200 300 400 500 600 700 800 900 

Time, (seconds) 


Figure 16. Tracker Comparison With One Sensor Type For The Non-Maneuvering 

Aircraft. 
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Test 1: Set (e) 

Tracking Algorithm Validity 

Description: 

Observe the performance of the trackers using the same sensor type with each of the trackers 

Sensor 

Airborne Surveillance 
Radar 

Identification Friend or 

Foe System 

Electronic Surveillance 
System 

Infrared Search 
and Track 
System 

Tracker 

PDAF 

IMM 

Kalman Filter 

CGKF 






Settings: 

1.0 

-SLMq^= 1.0 

- q"' = 1000.0 

-q^= 1000.0 


- 1 clutter pt 

-Turnq-= 10000.0 



Targets 

Maneuvering Aircraft 

Non-Maneuvering Aircraft 

Slow-Moving Ship 

Fixed Site 


Mean Error (-) vs. Time, 1 Sensor Type/4 Trackers/lOO run(s) 
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Figure 17. Tracker Comparison With One Sensor Type For The Slow-Moving Ship. 
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Test 1: Set (e) 

Tracking Algorithm Validity 

Description: 

Observe the performance of the trackers using the same sensor type with each of the trackers 

Sensor 

Airborne Surveillance 
Radar 

Identification Friend or 

Foe System 

Electronic Surveillance 
System 

Infrared Search 
and Track 
System 

Tracker 

PDAF 

IMM 

Kalman Filter 

CGKF 






Settings: 

1000.0 

-SLMq^= 1.0 

- q"' = 1000.0 

-q^= 1000.0 


- 1 clutter pt 

-Turnq’= 10000.0 



Targets 

Maneuvering Aircraft 

Non-Maneuvering Aircraft 

Slow-Moving Ship 

Fixed Site 


8000 

7000 
6000 |j 
5000 


Mean Error (-) vs. Time, 1 Sensor Type/4 Trackers/100 run(s) 


INFRARED SEARCH AND 
TRACK SIMULATOR 


PDAF 


Mean Measurement Error 

\ 



Figure 18. Tracker Comparison with One Sensor Type for the Fixed Site Target. 
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B. TEST 2: MULTIPLE SENSOR MULTIPLE TRACK 

For Test 2, the multiple sensor/multiple track case was examined to demonstrate 
the different capabilities of sensors and trackers. Representative plots for the 
maneuvering target were included. The sensors and trackers were paired up as shown in 


Table 5. 


Airborne Surveillance Radar 

Probabihstic Data Association Filter 
(PDAF) 

Identification Friend or Foe Detection 
System (IFF) 

Interacting Multiple Models (IMM) 

2-model filter 

Electronic Surveillance Measures System 
(ESM) 

Kalman Filter (KF) 

Infrared Search and Track System (IRST) 

Constant Gain Kalman Filter (CGKF) 


Table 5. Sensor/Tracker Pairs. 


Figure 19 is a display of the mean measurement errors in each of the sensors, 
while Figure 20 is the tracking solution for each sensor for the same run. Again, this test 
was simply a demonstration of the diversity in sensor capabihties. By observation, the 
trackers are sensitive to the amount of error in the sensors. For example, the IMM 
solution tracked very well due to the low errors in the IFF sensor. 

Some fine tuning in the filters could help improve tracking performance. 
However, a trade-off exists when adjusting the noise term, q . If the noise term is 
increased to reduce errors in the turn, the SLM performance will suffer as shown in Test 
1. This concept will be demonstrated again in Test 3, Set (c) by increasing the noise term 
in each tracker. 
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Test 2: Set (a) 

Multiple Target/Multiple Sensor 

Description: 

Demonstrate the differences between sensors, their associated mean measurement errors, and the resulting 
tracking solution. 






Sensor 

Airborne Surveillance 
Radar 

Identification Friend or 

Foe System 

Electronic Surveillance 
System 

Infrared Search 
and Track 
System 

Tracker 

PDAF 

IMM 

Kalman Filter 

CGKF 





Settings: 

1000.0 

-SLMq^= 1.0 

- q^= 1000.0 

-q^= 1000.0 


- 1 clutter pt 

-Turnq^ = 10000.0 



Targets 

Maneuvering Aircraft 

Non-Maneuvering Ahcraft 

Slow-Moving Ship 

Fixed Site 


Mean Measurement Error (-) vs. Time, 4 Sensors/100 run(s) 

6000 

5000 

4000 

3000 

2000 

1000 

0 

0 100 200 300 400 500 600 700 800 900 

Figure 19. Mean Measurement Error for Varying Sensor Types. 
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Test 2: Set (a) 

Multiple Target/Multiple Sensor 

Description: 

Demonstrate the differences between sensors, their associated mean measurement errors, and the resulting 
tracking solution. 






Sensor 

Airborne Surveillance 
Radar 

Identification Friend or 

Foe System 

Electronic Surveillance 
System 

Infrared Search 
and Track 
System 

Tracker 

PDAF 

IMM 

Kalman Filter 

CGKF 





Settings: 

1000.0 

-SLMq^= 1.0 

- q^= 1000.0 

-q^= 1000.0 


- 1 clutter pt 

-Turnq^ = 10000.0 



Targets 

Maneuvering Aircraft 

Non-Maneuvering Ahcraft 

Slow-Moving Ship 

Fixed Site 


Mean Track Error (-) vs. Time, 4 Sensors/4 Trackers/100 run(s) 



Figure 20. Tracking Solutions for the Multiple Sensor Case. 
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C. TEST 3: OBSERVATIONS ON ALGORITHM BEHAVOIR IN A SIMPLE 

COMBINED SOLUTION 

The purpose to Test 3 was to provide a demonstration of a multiple sensor system 
output, and to gain insight into the behavior of a combined tracker estimate. The 
combined solution was based on weights that were empirically derived from the results of 
Test 1, and the multiple sensor simulation results of Test 2. No additional data sources 
were used in determining the solution. The combined state result is simply the weighted 
average of each of the tracker solutions. This increases the amount of data stored by the 
program. For the, fusim simulation, the storage matrices and the function calls were a first 
step in implementing a track-to-track data association algorithm. Again, for this section, 
no additional filtering and no data association occurs in calculation of the weighted 
average of the track states. A true sensor fusion algorithm will be much more 
complicated [Ref. 1]. 

Combined Solution Set (a) Results 

The simple combined solution has initial conditions as described in Test 1 for the 
single sensor case. The output of the simulation is included in Figure 21. One sensor 
type is used with the four trackers to compare the trackers and to observe the results of 
combining the state outputs of each tracker in a rudimentary fashion. To further 
demonstrate the effects of diversity in the sensors and tracking algorithms, the Kalman 
filter was given a low noise term (q = 2) to allow an optimum solution in SLM. 

In the turn, only one of the 4 sensor/tracker pairs reported high errors (Figures 22 
and 23). This effect was by design due to the reduced noise factor, q in the Kalman 
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Filter. The overall error was significantly reduced by combining the results of the 
different systems. The non-maneuvering target (Figure 24) had the best results from the 
Kalman Filter, since the filter was optimized for straight line tracking. 

"IQ® Mean State (Tracking Solution) for All Sensors(1 run) 



LONGITUDE (feet) ^ .,^5 

Figure 21. Sample Plot of the Combined Solution. 


For the stationary target (Figure 25), the mean measurement error was greater due 
to the distance from the sensor platform. Again, the Kalman filter performed the best 
even though the sensor errors were high. The combined solution represents a 
compromise. 
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Test 3: Set (a) 

Simple Combined Solution for Multiple Target/Multiple Sensor 

Description: 

Demonstrate the results of combining the solutions into a single output state. 

Demonstrate the performance of the trackers using a single sensor type. 

Sensor 

Airborne Surveillance 
Radar 

Identification Friend or 

Foe System 

Electronic Surveillance 
System 

Infrared Search 
and Track 
System 

Tracker 

PDAF 

IMM 

Kalman Filter 

CGKF 





Settings: 

1000.0 

-SLMq^= 1.0 

- q^ = 2.0 

-q^= 1000.0 


- 1 clutter pt 

-Turnq^ = 10000.0 



Targets 

Maneuvering Aircraft 

Non-Maneuvering Ahcraft 

Slow-Moving Ship 

Fixed Site 


Mean Track Error (-) vs. Time, 4 Sensors/4 Trackers/100 run(s) 

6000 

5000 

4000 

3000 

2000 

1000 

0 

0 100 200 300 400 500 600 700 800 900 

Figure 22. Simple Combined Solution for Maneuvering Target. 
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Test 3: Set (a) 

Simple Combined Solution for Multiple Target/Multiple Sensor 

Description: 

Demonstrate the results of combining the solutions into a single output state. 

Demonstrate the performance of the trackers using a single sensor type. 

Sensor 

Airborne Surveillance 
Radar 

Identification Friend or 

Foe System 

Electronic Surveillance 
System 

Infrared Search 
and Track 
System 

Tracker 

PDAF 

IMM 

Kalman Filter 

CGKF 





Settings: 

q^ = 

1000.0 

- I clutter pt 

-SLMq*= 1.0 

^Turn^= 10000.|j 

- = 2.0 

-q^= 1000.0 

Targets 

Maneuvering Aircraft 

Non-Maneuvering Ahcraft 

Slow-Moving Ship 

Fixed Site 


Mean Track Error (-) vs. Time, 4 Sensors/4 Trackers/100 run(s) 

6000 

5000 

4000 

3000 

2000 

1000 

0 

0 100 200 300 400 500 600 700 800 900 



Figure 23. Simple Combined Solution for Maneuvering Target. 
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Test 3: Set (a) 

Simple Combined Solution for Multiple Target/Multiple Sensor 

Description: 

Demonstrate the results of combining the solutions into a single output state. 

Demonstrate the performance of the trackers using a single sensor type. 

Sensor 

Airborne Surveillance 
Radar 

Identification Friend or 

Foe System 

Electronic Surveillance 
System 

Infrared Search 
and Track 
System 

Tracker 

PDAF 

IMM 

Kalman Filter 

CGKF 





Settings: 

1000.0 

-SLMq^= 1.0 

o 

ri 

II 

-q^= 1000.0 


- 1 clutter pt 

- Turn = 10000.0 



Targets 

Maneuvering Aircraft 

Non-Maneuvering Aircraft 

Slow-Moving Ship 

Fixed Site 


Mean Track Error (-) vs. Time, 4 Sensors/4 Trackers/100 run(s) 



Figure 24. Simple Combined Solution for Non-Maneuvering Target. 
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Test 3: Set (a) 

Simple Combined Solution for Multiple Target/Multiple Sensor 

Description: 

Demonstrate the results of combining the solutions into a single output state. 

Demonstrate the performance of the trackers using a single sensor type. 

Sensor 

Airborne Surveillance 
Radar 

Identification Friend or 

Foe System 

Electronic Surveillance 
System 

Infrared Search 
and Track 
System 

Tracker 

PDAF 

IMM 

Kalman Filter 

CGKF 





Settings: 

1000.0 

-SLMq^= 1.0 

- q^ = 2.0 

-q^= 1000.0 


- 1 clutter pt 

-Turnq-= 10000.0 



Targets 

Maneuvering Aircraft 

Non-Maneuvering Aircraft 

Slow-Moving Ship 

Fixed Site 



Figure 25. Simple Combined Solution for Fixed Site. 
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Combined Solution Set (b) Results 


For the 2"‘^ set of runs, four different sensor types were used to observe the 
benefits of combining the data to form a composite tracking solution. Because of the 
diverse nature of the sensors, this implementation was somewhat representative of an 
actual sensor platform system. The combined solution was an attempt to take advantage 
of the best qualities of each of the sensor/tracker pairs. The model has a sensor that can 
handle clutter, a sensor that was designed to handle turning motion, an optimum straight- 
line-motion tracker, and a simphfied tracker. A simple truth that has been alluded to in 
previous sections is the nature of the Electronic Surveillance Measures (ESM) and the 
Infrared (IR) sensors. They are ESM and IR sensors in name only. In fact, the way that 
the ESM and IR sensors are modeled in the fusim program more closely resemble radar 
or lEE detection systems with greatly reduced resolution. Eor one other observation, it is 
unreahstic to expect an ESM system to maintain track on a moving/maneuvering target, 
since it requires a cooperative target. The ESM solution is for demonstration purposes 
only. 

Eor this set of runs (Eigures 26 through 29), the sensors did not have significantly 
large errors. However, the plots stiU show a significant reduction overall for the four 
sensors. In the turn, only one of the four sensors/tracker pairs reported high errors. 
Between the two maneuvering target Eigures (26 and 27), the differences in the plots 
show how target orientation with respect to the sensor platform can have an impact on the 
target errors and the combined solution. The error was significantly reduced by 
combining the results of the different systems. 
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As evidenced in all the previous plots (and previous tests), it would be 
advantageous to disregard the solution obtained with the CGKF until the initial tracking 
error has settled to below the measurement noise. This is another area where errors in the 
combined solution could be truncated. 

To demonstrate that the combined solution was not terribly susceptible to more 
reahstic errors in the IRST system, the errors were adjusted to observe the effects on the 
combined solution (Figures 30 through 31). 

In the resulting plots of Figures 30 and 31, the IRST/CGKF pair has an even 
larger initiahzation error due to track start-up. This unwanted result had an even greater 
effect on the combined solution, indicating a need to let the initiahzation errors settle 
prior to using the data with any confidence. 
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Test 3: Set(b) 

Simple Combined Solution for Multiple Target/Multiple Sensor 

Description: 

Demonstrate the results of combining the solutions into a single output state. 

Demonstrate the performance of the trackers using sensors of varying type 






Sensor 

Airborne Surveillance 
Radar 

Identification Friend or 

Foe System 

Electronic Surveillance 
System 

Infrared Search 
and Track 
System 

Tracker 

PDAF 

IMM 

Kalman Filter 

CGKF 





Settings: 

1000.0 

-SLMq^= 1.0 

- q^ = 2.0 

-q^= 1000.0 


- 1 clutter pt 

- Turn q^= 10000.0 



Targets 

Maneuvering Aircraft 

Non-Maneuvering Ahcraft 

Slow-Moving Ship 

Fixed Site 



Figure 26. Simple Combined Solution with Diverse Sensors, Maneuvering Target. 
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Test 3: Set{b) 

Simple Combined Solution for Multiple Target/Multiple Sensor 

Description: 

Demonstrate the results of combining the solutions into a single output state. 

Demonstrate the performance of the trackers using sensors of varying type 






Sensor 

Airborne Surveillance 
Radar 

Identification Friend or 

Foe System 

Electronic Surveillance 
System 

Infrared Search 
and Track 
System 

Tracker 

PDAF 

IMM 

Kalman Filter 

CGKF 





Settings: 

1000.0 

-SLMq^= 1.0 

- q^ = 2.0 

-q^= 1000.0 


- 1 clutter pt 

-Turnq-= 10000.0 



Targets 

Maneuvering Aircraft 

Non-Maneuvering Aircraft 

Slow-Moving Ship 

Fixed Site 



PDAw/ 
Surv Radar 


CGKF w/ 
IR sim \ 






Mean Track Error (-) vs. Time, 4 Sensors/4 Trackers/100 run(s) 
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0 100 200 300 400 500 600 700 800 900 


MANEUVERING TARGET 


Kalman w/ 
ESM Track 


Mean Meas Noise 
of Surv Radar 


Combined 

Solution 


Figure 27. Simple Combined Solution with Diverse Sensors, Maneuvering Target. 
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Test 3: Set{b) 

Simple Combined Solution for Multiple Target/Multiple Sensor 

Description: 

Demonstrate the results of combining the solutions into a single output state. 

Demonstrate the performance of the trackers using sensors of varying type 






Sensor 

Airborne Surveillance 
Radar 

Identification Friend or 

Foe System 

Electronic Surveillance 
System 

Infrared Search 
and Track 
System 

Tracker 

PDAF 

IMM 

Kalman Filter 

CGKF 





Settings: 

1000.0 

-SLMq^= 1.0 

o 

ri 

II 

-q^= 1000.0 


- 1 clutter pt 

- Turn = 10000.0 



Targets 

Maneuvering Aircraft 

Non-Maneuvering Aircraft 

Slow-Moving Ship 

Fixed Site 


Mean Track Error (-) vs. Time, 4 Sensors/4 Trackers/100 run(s) 



Figure 28. Simple Combined Solution with Diverse Sensors, Non-Maneuvering Target. 
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Test 3: Set{b) 

Simple Combined Solution for Multiple Target/Multiple Sensor 

Description: 

Demonstrate the results of combining the solutions into a single output state. 

Demonstrate the performance of the trackers using sensors of varying type 






Sensor 

Airborne Surveillance 
Radar 

Identification Friend or 

Foe System 

Electronic Surveillance 
System 

Infrared Search 
and Track 
System 

Tracker 

PDAF 

IMM 

Kalman Filter 

CGKF 





Settings: 

1000.0 

-SLMq^= 1.0 

o 

ri 

II 

-q^= 1000.0 


- 1 clutter pt 

-Turnq-= 10000.0 



Targets 

Maneuvering Aircraft 

Non-Maneuvering Aircraft 

Slow-Moving Ship 

Fixed Site 


6000 


Mean Track Error (-) vs. Time, 4 Sensors/4 Trackers/100 run(s) 
STATIONARY SITE 


5000 



Figure 29. Simple Combined Solution with Diverse Sensors, Fixed Site. 
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Test 3: Set{b) 

Simple Combined Solution for Multiple Target/Multiple Sensor 

Description: 

Demonstrate the results of combining the solutions into a single output state. 

Demonstrate the performance of the trackers using sensors of varying type 






Sensor 

Airborne Surveillance 
Radar 

Identification Friend or 

Foe System 

Electronic Surveillance 
System 

Infrared Search 
and Track 
System 

Tracker 

PDAF 

IMM 

Kalman Filter 

CGKF 





Settings: 

1000.0 

-SLMq^= 1.0 

o 

ri 

II 

-q^= 1000.0 


- 1 clutter pt 

-Turnq-= 10000.0 



Targets 

Maneuvering Aircraft 

Non-Maneuvering Aircraft 

Slow-Moving Ship 

Fixed Site 


Mean Track Error (-) vs. Time, 4 Sensors/4 Trackers/100 run(s) 



Figure 30. Combined Solution with Diverse Sensors, Altered Sensor Errors. 
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Test 3: Set{b) 

Simple Combined Solution for Multiple Target/Multiple Sensor 

Description: 

Demonstrate the results of combining the solutions into a single output state. 

Demonstrate the performance of the trackers using sensors of varying type 






Sensor 

Airborne Surveillance 
Radar 

Identification Friend or 

Foe System 

Electronic Surveillance 
System 

Infrared Search 
ancfTrack 
System 

Tracker 

PDAF 

IMM 

Kalman Filter 

CGKF 





Settings: 

1000.0 

-SLMq^= 1.0 

o 

ri 

II 

-q^= 1000.0 


- 1 clutter pt 

- Turn = 10000.0 



Targets 

Maneuvering Aircraft 

Non-Maneuvering Aircraft 

Slow-Moving Ship 

Fixed Site 


Mean Track Error (-) vs. Time, 4 Sensors/4 Trackers/100 run(s) 
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Figure 31. Combined Solution with Diverse Sensors, Altered Sensor Errors. 
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Combined Solution Set (c) Results 


2 

For the third set of runs in the Test 3 configuration, the noise term (q ) was 
increased in each of the tracking filters to allow better performance when tracking 
maneuvering targets. Also, observing the overall effect on the combined solution was an 
objective of this test set. The tracker parameters were altered as shown in Table 6. 


Tracker: 

PDAF 

IMM 

KF 

CGKF 

Previous q^: 

1000 

SLM: 1 

Turn: 10000 

2 

1000 

Set 3 q^: 

5000 

SLM: 1 

Turn: 100000 

5000 

5000 


Table 6. Altered Parameters of Combined Solution Set 3. 


For this set of runs, the overall combined solution was shghtly better than 
expected. Part of the reason for the success of the solution was the continued effective 
performance of the IMM. Although the IMM was paired with a more accurate sensor, 
increasing the turn q did not greatly affect the SLM performance. The Kalman filter 
solution for the ESM system showed the most dramatic increase in overall error (figure 
32), while the CGKF solution for the IRST system was only shghtly increased. 
However, the turn performance in each of the trackers improved, resulting in virtuaUy no 
increase in error due to turning motion in the combined solution. 

Against the non-maneuvering target (Figure 33), the ESM/Kahnan filter solution 

was no better than the IRST/constant gain case, even though the IRST sensor had much 

more error than the ESM system. Part of this poor performance in the KE can be 

attributed to a very large distance from the sensor platform. However, when compared to 

the low noise case (q^ = 2) of figure 31, the KE did a fine job of tracking over great 
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distances even when the sensor error was fairly high (as in the case of the ESM system). 
The decreased performance of the Kalman filter from Figure 31 to Figure 33 is directly 
attributable to the increased noise term. The Kalman filter also had poor performance 
against the fixed site in Figure 34. The poor performance of the KF is an indication that a 
compromise is necessary when using a standard Kalman filter for tracking. Should the 
Kalman filter be optimized for SFM performance, or for turning performance? The 
Kalman solution with a low noise term is the optimum solution when the target is not 
maneuvering. In Test 3a, the Kalman filter had the least error against the non¬ 
maneuvering target and the fixed site (a same sensor comparison). 
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Test 3: Set (c) 

Simple Combined Solution for Multiple Target/Multiple Sensor 

Description: 

Demonstrate the results of combining variable solutions into a single output state. 

Demonstrate the performance of the trackers with increased noise terms. 






Sensor 

Airborne Surveillance 
Radar 

Identification Friend or 

Foe System 

Electronic Surveillance 
System 

Infrared Search 
and Track 
System 

Tracker 

PDAF 

IMM 

Kalman Filter 

CGKF 





Settings: 

- = 5000.0 

-SLMq*= 1.0 

- q^ = 5000.0 

- q^ = 5000.0 


- 1 clutter pt 

-Turnq^= 100000.0 



Targets 

Maneuvering Aircraft 

Non-Maneuvering Aircraft 

Slow-Moving Ship 

Fixed Site 


CD 

(D 


O 

i5 


6000 


5000 


+ 

■Nf- 


Mean Track Error (-) vs. Time, 4 Sensors/4 Trackers/100 run(s) 

MANEUVERING 

TARGET 



Time (seconds) 


Figure 32. Combined Solution with Diverse Sensors, Altered Tracker Parameters. 
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Test 3: Set (c) 

Simple Combined Solution for Multiple Target/Multiple Sensor 

Description: 

Demonstrate the results of combining variable solutions into a single output state. 

Demonstrate the performance of the trackers with increased noise terms. 






Sensor 

Airborne Surveillance 
Radar 

Identification Friend or 

Foe System 

Electronic Surveillance 
System 

Infrared Search 
and Track 
System 

Tracker 

PDAF 

IMM 

Kalman Filter 

CGKF 





Settings: 

- = 5000.0 

-SLMq^= 1.0 

- q’ = 5000.0 

- q^ = 5000.0 


- 1 clutter pt 

- Turn = 100000.0 



Targets 

Maneuvering Aircraft 

Non-Maneuvering Aircraft 

Slow-Moving Ship 

Fixed Site 


Mean Track Error (-) vs. Time, 4 Sensors/4 Trackers/100 run(s) 



Figure 33. Combined Solution with Diverse Sensors, Altered Tracker Parameters. 
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Test 3: Set (c) 

Simple Combined Solution for Multiple Target/Multiple Sensor 

Description: 

Demonstrate the results of combining variable solutions into a single output state. 

Demonstrate the performance of the trackers with increased noise terms. 






Sensor 

Airborne Surveillance 
Radar 

Identification Friend or 

Foe System 

Electronic Surveillance 
System 

Infrared Search 
and Track 
System 

Tracker 

PDAF 

IMM 

Kalman Filter 

CGKF 





Settings: 

- = 5000.0 

-SLMq^= 1.0 

- q’ = 5000.0 

- q^ = 5000.0 


- 1 clutter pt 

-Turnq-= 100000.0 



Targets 

Maneuvering Aircraft 

Non-Maneuvering Aircraft 

Slow-Moving Ship 

Fixed Site 


6000 


5000 
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a> 

a> 


3000^ 
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■+ 
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Mean Track Error (-) vs. Time, 4 Sensors/4 Trackers/100 run(s) 
FIXED SITE 


Mean Measurement 
Error for Surv Radar 
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LU 


Surv Radar 


Combined 

Solution 
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Figure 34. Combined Solution with Diverse Sensors, Altered Tracker Parameters. 


107 
























THIS PAGE INTENTIONALLY LEET BLANK 


108 



VI. CONCLUSIONS AND RECOMMENDATIONS 


Testing algorithms in a simulation environment can be an effective way to study 
the sensor fusion problem without using hve sensors and targets. The cost of trying to 
use actual targets with repeated initial conditions would be astronomical. A simulation 
scheme can test and retest an algorithm using randomly generated perturbations to 
provide an indication of the algorithm effectiveness. The, fusim program was written to 
give the user flexibihty in selecting a sensor platform, the sensors used on the platform, 
the target types, the problem orientation, and the tracking algorithms to be used with the 
sensors. To change parameters, the user can simply alter the MATLAB code within the 
applicable files to produce the desired output. The fusim program can be used to compare 
tracking algorithms in a multiple sensor/multiple target environment, view the effects of 
using diverse sensors, and evaluate the effects of data association algorithms. For this 
early edition of the fusim program, only a rudimentary combination of sensor/tracker 
outputs was used to evaluate the effects of multiple, diverse sensors. 

The fusim program was designed to allow averaging over numerous runs to find a 
mean solution for each sensor and track update. The error generation was based on 
Gaussian distributed random white noise, meaning that numerous runs and average 
solutions would be necessary to account for the randomness of the noise. For most test 
cases, 100 or 200 runs were used to produce a smooth track output and error plot. 
However, this information tells nothing of the tracker performance for a single run. 
Therefore, the program also allows plotting results from a single run. Data gathered over 
a single run may be more relevent for evaluating a tracking or fusion algorithm once the 
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algorithm is proven. When used against hve targets, the tracker wiU have to perform 
consistently and predictably. 

In this study, the effect of multiple sensors in different locations contributing to 
the tracking picture was not analyzed. None of the track results were able to take 
advantage of multiple sensor locations, as in the case of a Cooperative Engagement 
Capabihty system with multiple participating units. Also, none of the tracking algorithms 
were true multiple sensor/multiple track algorithms. Throughout the simulation, every 
measurement was assumed to be associated with only one track. To create a true multiple 
sensor tracking routine would require measurement-to-track and track-to-track data 
association. Handhng of multiple tracks would require the use of a Multiple Hypothesis 
Tracker or a Joint Probabihstic Data Association Filter, or some other means of 
measurement-to-track data association. 

The performance of each of the trackers was highly dependent on the errors 
within each sensor, and in the way that each sensor measured target data. One of the 
artificial aspects of the fusim program is that all sensors report measurements in range 
and bearing. For example, the modeling of the ESM and IR sensors was greatly 
oversimplified in the fusim program to allow simple application of the measurement and 
tracking functions. Instead of using bearing only, such as in the case of an ESM target, 
error was added to an x-y position report. For all sensors in the fusim program, errors 
were modeled in range and bearing. The true x-y position reports were converted to 
range and bearing, the range/bearing errors were added, the position report was converted 
back to x-y coordinates, and the ‘noisy’ measurement (and error covarinance) was 
returned to the main program for use in the tracking functions. This same procedure 
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occurred for all sensors and tracks regardless of type. In Test 1 Set (a), where the sensor 
errors were zero, the tracker that was used was irrelevant for aU cases except the 
maneuvering target. Some error was noted in the turn, suggesting that the kinematics of 
the problem has an effect on the target state even if the errors in the sensor are zero. 
Also, nearly every test case involved a moving sensor platform. The sensor platform 
implementation was a little more reahstic, but the analysis was more challenging since 
the orientation of the problem changed continuously. The data is not as smooth as it 
would be for a fixed sensor platform with no navigation errors in the position. This 
analysis provided a general idea of how each of the trackers performed from the 
perspective of the sensor platform. 

The fusim program was used to compare the performance of several different 
tracking algorithms. The PDAF was the only tracker in this evaluation that was forced to 
deal with clutter. Because of the clutter variable, the PDAF was not as consistent and 
predictable as the other trackers. On one set of runs, the PDAF ran away from the actual 
position of the stationary target due to a shght velocity induced in the state vector by the 
clutter point. Since the clutter point was placed around the predicted state and the tracker 
does not distinguish between the clutter measurement and the actual measurement, the 
state was allowed to move away from the actual position of the target. The data link 
report followed suit, and provided a good illustration of how easy it is to flood a system 
with bad information. The information output from a system is only as good as the 
information input. 

The IMM was the most consistent tracker of the evaluation set for maneuvering 

targets. However, the IMM also performed fairly well with non-maneuvering targets, 
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slow-moving targets, and fixed sites. In the algorithm, the only difference between the 
two models was the noise factor q in the SLM model and the turn model. The IMM 
produced the best overall tracking solution every time, but required much more code and 
computation time (see Appendix B). 

The Kalman filter performed very well when tracking non-maneuvering, slow- 
moving, and fixed-site targets. The Kalman filter did not handle maneuvering targets 
very well. By increasing the plant noise, the Kalman filter handled the turns better, but 
resulted in noisier straight-hne tracking with a higher average error. 

Repeatedly, the Constant Gain Kalman Filter produced reasonable straight-line 
performance, but did not do as well during initiahzation or in turning motion. The 
tracker required the fewest computations, but the large track initiation errors were present 
in every test case, except for the zero-error case. However, when the measurement errors 
were low due to sensor accuracy or due to the target position relative to the sensor 
platform, the CGKF performed as well as the PDAF or the IMM for the same 
sensor/target pair. This suggests that the CGKF would be useful in high-workload 
situations where non-priority tracks could be tracked with less computational strain, or 
the CGKF would be useful when sensor errors are low. 

The applications of sensor fusion are widely varied, but particular attention is paid 
to the military apphcation in this thesis. Battle group operations, the Common 
Operational Picture, land attack scenarios, time-critical strike, and sateUite data 
apphcations are just a few examples of the military uses of Sensor Fusion. The ‘Sensor 
Fusion’ problem is very broad in scope, to include network centric interconnectivity. 
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platform level processing and display, and global adaptivity for flexible mission 
application. This thesis does not address the global aspects of sensor fusion such as 
interoperabihty or systems engineering, nor does it address platform level systems 
engineering. Rather, this thesis was an attempt to reach into the core of sensor fusion 
algorithms and discuss the tracking and data association aspects of tracking algorithms. 

The fusim program was designed with a great deal of flexibihty to allow for future 
apphcations and upgrades. It is recommended that this study be continued to allow the 
following expansions: 

1. Extend this 2-dimensional study to 3-D with full coordinate system 
accounting. 

2. Produce functions that allow full data association for (1) measurement-to- 
track (2) track-to-track (3) data-to-track. 

3. Integrate an emitter library, similar to the sensor file, that could be used to 
provide more reahstic ESM system inputs, and attribute information. 

4. Expand the program to include full accounting of track attributes for 
determining track identification, and allow attribute tracking. 

5. Design the Graphical User Interface for max im um user benefit of fusim 
and study the Human-Machine Interface aspects of Sensor Eusion. 
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APPENDIX A. FUSIM FUNCTIONAL DESCRIPTION 


Fusim Functional Description 
fusim 

The FUSIM file is the main control file for the entire simulation. FUSIM calls the 
various associated functions to set up the needed matrices and variables, creates the truth 
data for the targets, adds the measurement error, sends the noisy measurements to the 
tracker, and then stores and plots the target states and associated least square errors. The 
FUSIM function is called from the user interface or from the MATLAB command hne. 
The MATLAB code for aU the functions is included in Appendix B. First, the START 
and FTURN functions are called to initialize aU the timing and target parameters, as input 
by the user. From this, the total number of maneuvering and non-maneuvering targets are 
determined and the sensor information (including number and types of sensors) is 
extracted. Based on the user selection, the primary tracker type is extracted for the 
primary sensor, and the other sensors are assigned the remaining trackers. This 
information is stored in the (trkr) vector. 

Next, the outer loop for the Monte Carlo runs is set up, and all appUcable 

variables and matrices become initialized. This way, for repeated runs, the variables are 

reset to nuU or to their initial values at each iteration of the loop. The generation of truth 

data begins with extracting the initial positions using the NONMAN and MANMATRIX 

functions. The maneuvering and non-maneuvering targets are kept separate for 

generating the truth data, and then regrouped for generating the noisy measurements. 

The initial position and velocity vectors look l ik e [xl;vxl;yl;vyl, x2;vx2;y2;vy2, .] 
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for each target. In order to generate the second measurement, the TIMESTEP function is 
called to iterate the target motion in time. 

The MANMATRIX function has an additional task that is not part of the 
NONMAN function. MANMATRIX also returns the (mantime) matrix that contains the 
column vectors of maneuver times for the maneuvering targets only. However, for the 
first two measurements, straight hne motion for all targets is assumed. The truth data is 
then stored in the (posout) matrix as x-y pairs stacked for each target. A new column of 
x-y pairs for each target is created with each time step. 

The first two measurements are required for track initiahzation. This is a major 
assumption of this simulation. AU tracks are initiahzed with a Eeast Squares Estimation 
track initialization, regardless of the tracker used in the simulation. Trackers that use 
Kalman predictions and updates require an initial state and covariance for the tracker to 
work. The INIT function is called for each sensor of the sensor platform. This is where 
the multi-dimensional capability of MATEAB is first exercised. While the 2-D form of 
the matrix is needed for the calculations, storage of the multiple sensor data is 
accomplished by adding a “page” to the 2-D matrix for each additional sensor. This 
gives the matrices depth, or a third dimension. Storing this data separately is necessary 
because each sensor has a unique accuracy. The SETZTRUE function is used to 
conveniently stack the z-true data for storage and plotting. The noisy measurements are 
also saved for error comparisons after the tracker has updated the state (zout). 

At this point, the simulation is completely initialized and ready to start looping 
through the main body of the program. Eor each iteration, the following actions occur: 
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1. Truth data is generated for non-maneuvering targets with TIMESTEP. 


2. A test is performed for each maneuvering target to see if the target is in 
straight hne motion or in a maneuvering state. If the target is in SEM, the straight line 
(E) matrix is sent to TIMESTEP. If the target is in a tuning state (based on ‘mantime’ 
matrix times), the (Eturn) matrix is sent to TIMESTEP for proper modehng of the turn. 

3. AU of the truth data is combined and sent to the MEAS function for generation 
of noisy measurement data. The error and the error covariance is returned for each target 
and for each sensor, which add to the 3-D matrices already generated by INIT. 

4. The noisy measurement data, error covariance, initial state and initial state 
covariance are aU sent to the applicable tracker for each sensor. This data is aU kept 
separate for plotting the individual solutions and errors. 

5. Once aU the data and tracking solutions have been calculated for the entire 
simulation run, the ESECAEC function determines the Eeast Square Error between the 
tracker state and the truth position of each point and for each sensor. If multiple Monte 
Carlo runs are conducted, the results are stored and averaged over the number of runs. 
This averaged data is used for plotting the final results. An additional option is allowed 
in the code for plotting the state solutions at each loop iteration (SIMPEOT). The 
functions TRACKPEOT and ERRORPEOTS are used to plot the results. 

When running/w5/m, the track plot will look similar to the results in Eigure 35. 
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X 10* Mean State (Tracking Solution) for All Sensors(20 runs) 



LONGITUDE - in feet for now ^ 


Figure 35. Representative PPL 


[A,Al,delta,tjnloopSjnsampleSjSpi,flag,ownsens] = start(dumb); 

The START function is called first by FUSIM to set up the initial variables 
needed for the simulation. The START function reads the START.MAT file to collect 
timing information about the simulation, initiates the matrices of targets, and initiates 
other needed matrices (F,H). The notable functions and variables of the start function are 
described in the following sections. 

The START function begins by reading the START.MAT file to initialize the 
timing variables and other user-defined variables. The START.MAT information 
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initializes variables of the (simstart) vector in the following format: [t(simtkne), nloops, 
deltatkne, DLRP, placeholder, tracker selection 1-4]. The total simulation time (t) is 
input in minutes and then converted to seconds. The number of loops (nloops) 
determines the total number of Monte Carlo runs to be executed by the program. A large 
number of runs can add significantly to the execution time required by the simulation. 

The time variable (deltatime) is the update rate of the sensors given in tenths of 
minutes. For example, if 0.1 is used, 0.1*60 = 6 seconds. In other words, the update rate 
of the primary sensor is every 6 seconds. The variable (DLRP) is used to center the 
entire simulation at some location on the earth. When entered in the simulation, (DLRP) 
becomes the origin of the X-Y coordinate grid. Note: The (DLRP) variable is not used 
in this version of fusim. The tracker selection is used in a later function. Now that the 
total simulation time and the time step (delta) are known, the total number of samples 
(nsamples) can be calculated. The variable (nsamples) is an integer value that provides 
FUSIM with the total number of samples or measurements to be taken for each of the 
targets. 

Next, the TARDAT.MAT file is read into the (Al) matrix. By examining the 
elements of the (Al) matrix, a determination is made whether the target is maneuvering 
or not. This information will be used later for setting up the truth data. As the maneuver 
test is performed, the targets are separated into the (A) matrix (maneuvering targets) and 
the (Al) matrix (non-maneuvering targets). These two matrices contain all the truth data 
about all the targets as input by the user. The truth elements of the sensor platform and of 
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the target vectors are nearly the same, as shown below and further described in Appendix 


B: 


Example sensor platform vector (format of TARDAT.MAT data): 

[N,365600,W,0761700,110,175,25000,A/C,e2c,radarl,iff,esm,ir,comm,legtime,L/ 
R,turn time] 

1: North or South Latitude N = 0, S = 1 (not utihzed in this version) 

2: X position, with 0/0 as the origin (feet) 

3: East or West Longitude: W = 0, E = 1 (not utihzed in this version) 

4: Y position, with 0/0 as the origin (feet) 

5: Course in degrees 
6: Speed in Knots 

7: Altitude in feet (not utihzed in this version) 

8: Type: 1 = Aircraft, 2 = ship, 3 = ? (not utihzed in this version) 

9: Name 1 = E-2C, 2 = ? (not utihzed in this version) 

10: Sensor 1 (sensor platform) / Emitter 1 (target) 

11: Sensor 2 (sensor platform) / Emitter 2 (target) 

12: Sensor 3 (sensor platform) / Emitter 3 (target) 

13: Sensor 4 (sensor platform) / Emitter 4 (target) 

14: communications type (not utihzed in this version) 

15: Leg Time - how long the legs of the orbit will be (if turning) 

16: Turns - Left or Right, 0 = Left Turn, 1 = Right Turn 
17: Turn Time - how long the turns wih be (if turning) 

Three other functions are initiahzed for the sensor platform. The sensor platform 
position (Spi) is initiahzed for use later in the program, flag is set to 1 if the sensor 
platform is maneuvering, and the vector of sensor types (ownsens) is initiahzed with the 
user-defined selections of sensor types. Up to 4 different sensors can be selected for the 
sensor platform. 
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Two other matrices are initialized in the START function, the (H) matrix and the 


(F) matrix. These are declared as global variables because they are used by several other 
functions in the FUSIM program. The (H) matrix is used simply to extract the x-y 
position variables from the state vector. The (F) matrix is used to increment in time a 
state vector for predictions and for estabhshing truth data when targets are in straight-line 
motion. 


[Fturn,B] = fturn(A,delta); 

The function FTURN is called next by the FUSIM control function to set up the 
needed matrices for maneuvering targets. An fturn matrix is used by each maneuvering 
target to model the target motion when the target is in a turning state. The fuU (Fturn) 
matrix is returned with the stacked 4x4 (Fturn) sub-matrices for the maneuvering targets 
only. [Fturn(ownship);Fturn(targetl);Fturn(target2);....]. The (B) matrix is set up by 
determining whether the target is turning left or right, extracting and converting the 
heading and speed information, and the initial truth state vector for each maneuvering 
target. The data for target 1 (usually ownship) occupies column 1, target 2 is in column 2 
and so on for the total number of maneuvering targets. 


function trkr = setuptrackers(ns) 

The function SETUPTRACKER uses the input value (ns) to set up a vector of 
sensor types. The value (ns) is based on the number of sensor types selected by the user. 
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The only tracker that the user can select is the primary tracking function desired for the 
primary sensor. The rest of the sensors are arbitrarily assigned a tracking function. 
Currently, the trackers are assigned as follows: 

1 = Probabilistic Data Association Filter 

2 = Interacting multiple Models 

3 = Kalman 

4 = Const Gain Kalman 

function [mantime,xim] = manmatrix(A,B,t) 

The function MANMATRIX constructs the Maneuver Time Matrix based on the 
size of the largest maneuver vector. First, the function extracts the maneuver times from 
the (A) matrix and converts the user-defined times to seconds. Using the leg and turn 
times, the maneuver vectors are then created with the MANEUVERTIME function, as 
discussed in the next section. The maneuver times for target I occupies the first column, 
the times for target 2 occupy the second column, and so on. If the vectors are different in 
length, the number of rows in the maneuver matrix is based on the longest vector. The 
rest of the matrix is filled in with zeros. The MANMATRIX function also returns the 
initial truth position (xim) for the maneuvering targets. 
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function [mt] = maneuvertime(t,legt,turnt) 


The function MANEUVERTIME is called for each target that is maneuvering. 
Based on the total time of the simulation and the leg and turn times of the target, the 
function creates a vector of cumulative maneuver times. The first time element in the 
vector wiU always be the first leg time. The next time will be the length of time required 
to complete the first turn added to the first leg time. The target uses the straight-line 
motion model (E) for the first leg, and the turning model (Eturn) during the time that the 
target is turning. After completion of the turn, the target returns to SEM until the next 
turn time is reached. As an example, if the first leg is 600 seconds, the target will use 
SEM model for 600 seconds. If the turn takes 30 seconds to complete, the turn model 
wiU be used from time 601 to time 630. At time 631, the target returns to SEM until time 
631 - 1 - 600 or 1231 seconds. This target motion model continues until the total time of the 
simulation is reached. 


function [xin] = nonman(Al) 

The non-maneuvering targets are initiahzed using the NONMAN function. This 
function extracts the position and velocity data from the (Al) matrix of non-maneuvering 
targets, converts the speed to feet per second, and breaks down the velocity into x and y 
components based on the heading. This also includes fixed targets that have no motion 
associated with them. The state vector is constructed for each target and stored as a 
column of the non-maneuvering matrix. 
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function [spxi] = ownxi(ownship) 


This function simply extracts the initial position and velocity data from the 
ownship initial parameters vector (as input by the user). 


function [Q,P,R,derror,xhat,zret] = init(allxl,allx2,delta,Spi,Sp,sensor) 

The INIT function is one of the most important functions in the simulation. This 
function performs nearly all the same operations as the main FUSIM loop. Since all the 
trackers used in this simulation require some form of initialization, a simple Least 
Squares Estimation initialization is performed for each. Using the least squares 
initialization for all the trackers is a bit of a simplification, but the simulation is much 
more flexible. Nearly any tracker that uses a Kalman type prediction and update can be 
plugged into this simulation and evaluated. The INIT function uses the first two noisy 
measurements (assuming SLM) to initialize the Kalman filter state and the prediction 
covariance. Using the initial truth data provided by the user, the function POLE2CART 
is called with the apphcable sensor range and bearing error, and the error covariance 
matrix. Returned from POLE2CART are the noisy measurement, the measurement error, 
and the new error covariance (based on the measurement). More information on 
POLE2CART is provided in the following section. 

Since two measurements are required for the initiahzation, the targets are stepped 
through time for one time increment using the TIMESTEP function. The initiahzation 
assumes SLM for ah targets. The truth data is stored for later output, and POLE2CART 
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is called again to create a noisy measurement for each target based on the sensor errors. 
With the 2 noisy measurements, sufficient data exists for the least squares initiahzation. 
The data that are returned include the noisy measurements (zret), the initial states (xhat), 
the error covariance (R), the prediction covariance (P), the measurement errors (derror), 
and the plant noise (Q) matrix. 


function [z,R,disterror] = pole2cart(ztrue,Sp,sr,sb,Sv) 

The POLE2CART function simulates the errors inherent within a sensor by 
adding error to the truth position and returning a new noisy measurement. First, the x-y 
truth data is converted to polar coordinates based on the sensor position and the target 
position. Next, the range and bearing errors are multiphed by a randomly generated 
number to simulate Gaussian distributed white noise, and then added to the true range 
and bearing. The new range and bearing are then converted back to Cartesian 
coordinates and the error covariance matrix is constructed. The error between the noisy 
measurement and the truth position is measured and returned with the noisy position and 
error covariance matrix. Using POLE2CART for the measurements creates a shght bias 
error that is ignored for this simulation [Ref. 1]. 
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function [spP,spR,sperror,spstate,spmeas] = spinit(spxi,spx2,Q,D) 


The SPINIT function provides the initial gain and state estimation for ownship 
navigation, assuming straight line motion. The initialization is precisely the same as the 
target track initializations. 


function ztrue = setztrue(m,n,xin,xim) 

The SETZTRUE function extracts the true x-y data from the truth state vector. 
This is performed for both maneuvering and non-maneuvering targets. The data is stored 
in the (ztrue) vector and returned to the EUSIM function for storage in the posout matrix. 


function [sigr,sigb,Sv] = senserror(sensor) 

The SENSERROR function matches the type of sensor defined by the user for the 
sensor platform to the corresponding errors. Eor each sensor on the sensor platform, the 
range and bearing errors are returned along with the sensor variance matrix. Eor 
simphcity, all sensors are treated the same. The errors are given in range and bearing, 
allowing use of the POEE2CART function for all sensor measurements. This is not the 
most accurate representation of sensors, however, treating the sensors the same allows 
more flexibihty in the simulation. The following hst describes the sensor type and the 
errors used for the simulation: 
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Case 

Sensor Type 

? range, ft 

bearing) 

Sensor 

Variance 

1 

Airborne Surveillance Radar 

100 

.005 

99 2 0 

9 ■ 9 

^ ^ hearing f 

2 

Surface Surveillance Radar 

20 

.0001 

same as above 

3 

Airborne Fighter Radar 

20 

.001 

a 

4 

Surface Fire Control Radar 

10 

.0001 

a 

5 

Non - Cooperative Target 

Recognition 

20 

.001 

a 

6 

Eastern Identification Friend or Foe 

50 

.01 

a 

7 

IFF Interrogate 

50 

.001 

a 

8 

IFF Reply 

20 

.001 

a 

9 

Airborne Electronic Surveillance 
Measures (ESM) 

500 

.01 

a 

10 

Surface ESM 

300 

.01 

a 

11 

Infrared Detector 

6000 

.0001 

a 

12 

Ownship Navigation 

10 

.00001 

a 

13 

Data Link Tracking Errors 

200 

.01 

a 


Table A-1. Included Sensor Types and Associated Errors. 


function [xnew] = timestep(F,x) 

The TIMESTEP function receives the time increment matrix and returns the new 
state vector for a moving target. If the target is in SEM, the (E) matrix is used. If the 
target is turning, the (Eturn) matrix is sent. The TIMESTEP function simply multiplies 
the step matrix and the state vector. 


[spxtrue,spmeas,sper,spR,in]=spitmeas(sensor,spx2,spman,ownfturn,simt,in) 

The SPITMEAS function is used to iterate the sensor platform position (Sp) and 
to calculate the noisy measurement. Eirst, a test is conducted to determine if the sensor 
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platform is a maneuvering condition. If maneuvering, the (ownfturn) matrix is used to 


update the truth position of the platform. Assuming that the origin of the navigation 
solution is ownship eg, the previous value of Sp is used as the origin for the call to 
POLE2CART. 


function [R,zret,derror] = meas(sensor,zre,Sp) 

The MEAS function receives the x-y truth data as well as the sensor type and 
sensor position. A function call is made to the SENSERROR function to get the sensor 
errors for the current sensor. These errors are sent to the POEE2CART function for 
calculation of the noisy measurements. The function then stacks the noisy measurements 
in a vector (zret) and stacks the measurement covariance matrices (R) for return to the 
main program. This process is repeated for all targets and for each sensor type. 


function [xhat33,count] = tracker(Ql,PI,Rl,xhatl,zl,count,k,trkr) 

A different tracker is called for every sensor that has been selected for the sensor 
platform. The TRACKER function is sent all the apphcable covariance data, the 
previous state, the measurement, and the tracker type. To call the appropriate tracker, a 
switch-case statement is used with identical argument hsts for all available trackers. The 
tracker types are hsted with the SETUPTRACKER function. 
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function [xhat33,count] = pdaf(Q,P,R,xhat,zret,count) 


function [xhat33,count] = immtrk(Q,P,R,xhat,zret,count) 

function [xhat33,count] = kalman(Q,P,R,xhat,zret,count) 

function [xhat33,count] = constGKF(Q,P,R,xhat,zret,count) 

The PDAF function is the Probabilistic Data Association Filter, which can be 
used where clutter is present in the target tracking gate. The IMMTRK function contains 
the code for the Interacting Multiple Models tracker. All the needed variables are defined 
and initiahzed the first time the function is called. The variables that need to be saved 
from iteration to iteration are simply declared as global variables within the tracking 
function. This way, the variables are not destroyed each time the function returns the 
updated state. The MATLAB code is included in Appendix B. The function calls for the 
Kalman filter and the constant gain Kalman filter are identical. 


function [linkdat,lstate] = linkreport(xhat,simt) 

This function creates a simulated data hnk report for the targets and allows for 
addition of attribute data, such as Identification Friend or Foe data. The format of the 
report is very flexible and can be changed as desired by user to suit individual needs. 
[x,y,velocity,heading,time,iffl,iff2,iff3,iff4,iffc,emitter l,emitter2,ID,reportingunit] 

Within the function, random noise is added to the to velocity and position to 
simulate navigation errors between platforms. Also, to simulate DLRP/gridlock errors, a 
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fixed value of bias error is added to the position. The hnk report is updated every 5 delta¬ 


time cycles. 


function simplot(mn,ns,xhat) 

The SIMPLOT function allows plotting of the track data as the program is 
executing. If multiple Monte Carlo runs are used in execution of the simulation, this 
function is not recommended since execution is slowed dramatically. 


function LSE = lsecalc(mn,state,posout,nsamples) 

The function LSECALC calculates the Least Square Error between the state 
vectors and the truth positions. This data is stored for each target and each sensor at 
every time step. 


function handl = errorplots(mn,ns,nsamples,delta,merror,LSmean) 

The ERRORPLOTS function plots aU the measurement error and Least Squares 
error for every target and sensor vs. time. 
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function hand = trackplot(mn,ns,posout,zoutmean,statemean) 


The TRACKPLOT function plots the true position, the noisy measurements, and 
the track states for each track and each sensor. Also provided is a single plot for 
comparison of all sensors for each of the tracks. 

Running the Fusim Simulation 

To run fusim: Type fusim at the command line and use the default simulation 
time, sensor types, and targets. To change the simulation time, open the START.MAT 
file and change the first element to the desired simulation time (in minutes). The second 
element is the number of runs. The number of simulations should be set to 1 if using the 
SIMPLOT command in FUSIM. If numerous runs are desired for smoothing of the 
randomly generated error data, set the element to the desired number (recognizing that 
simulation time increases dramatically). To change the timestep, delta, set the third 
element to the desired delta in minutes. 

To change the primary tracker from PDA to IMM or Kalman or Constant Gain 
Kalman, change the final element (1) to a (2), (3), or (4) respectively. The assignment of 
trackers to other sensors will fall in the order given in the SETUPTRACKERS function. 
To change the sensor types in the sensor platform, open the TARDAT.MAT file and 
change elements 10-13 to the desired sensor types as described in the SENSERROR 
function. Einally, any number of targets and target types can be entered in the same 
format as described previously. 
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APPENDIX B. FUSIM PROGRAM MATLAB CODE 


FUSIM PROGRAM MATLAB CODE 


NAME 

Fusim 

INPUT 

None 

OUTPUT 

None 

PURPOSE 

Fhe. fusim program is the main simulation file. This file initializes aU variables, creates 
the simulation truth data, calls the measurement function, calls the trackers, calls the data 
link simulator, and calls the plotting functions. 


clear aU; 
format long; 


global F H Sp; 
dumb = 0; 

[A,A1,delta,t,nloops,nsamples,Spi,flag,ownsens] = start(dumb); 

%Determine if ownship platform (the sensor platform) is maneuvering and set up for 
% separate tracking 
spxi= []; 
if flag == 1 

% Ownship is maneuvering 
ownship = A(l,;); 

A(l,:) = []; 

[ownfturn,ownB] = fturn(ownship,delta); 

[spman,spxi] = manmatrix(ownship,ownB,t); 
elseif flag == 0 

% Ownship is not maneuvering 
ownship = Al(l,:); 

Al(l,:) = []; 

spxi = nonman(ownship); 

end 

[Fturn,B] = fturn(A,delta); 

% Fturn comes back with stacked Fturn matrices [Ftl;Ft2;Ft3;....] 
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% B comes back with all the other data stacked; [turng;hdg;speed;xi turng;hdg;speed;xi,...] 

% The total number of targets wiU always be n + m, or mn 
[n,b] = size(Al); 

[m,b] = size(A); 
mn = m+n; 
nsens = []; 

fork= l:max(size(ownsens)) 
if (ownsens(k) >= 1) 

nsens = [nsens,ownsens(k)]; 
end 
end 

ns = max(size(nsens)); % This gives the total number of sensors on the sensor platform, 
trkr = setuptrackers(ns); % Returns a vector with the desired trackers identified for %tracker.m 
% Try getting an average for a number of runs: 
for kk = Tnloops 

Sp = Spi; % Initializes the sensor position 

simt = 0.0; % Simt keeps track of the simulation time for the maneuvering targets 
posout = []; sppos = []; 
ztrue = []; 

zout = []; spzout = []; 
state = []; spstate = []; 
error = []; sperror = []; 
xin= []; 

x2n = []; spx2 = []; 
xim = []; 
x2m = []; 
xy=[]; 
xhatl = []; 

Ri = []; 
zi = []; 

Qi = []; 

PI = []; 

derrorl = []; 
count = 1; 
hnkout = []; 


138 



LSerror = []; 

% HRST TARGET TIMESTEP EOR INITIALIZATION: 

% Call the function 'nonman' to get the non-maneuvering target data in the same format 

% as the maneuvering target: xi will be a matrix of [x;vx;y;vy;'s.] side by side 

% mantime contains the column vectors of maneuver times for the maneuvering targets %onIy. 

sppos = [sppos,H*spxi]; 
spx2 = E*spxi; 
sppos = [sppos,H*spx2]; 
if n > 0 

xin = nonman(AI); 
x2n = zeros(size(xin)); 
for k = I :n 

x2n(:,k) = timestep(E,xin(:,k)); 
ztrue = [ztrue;H*xin(:,k)]; 
end 
end 

% Get the xi's, which wiU be used for truth data [x;vx;y;vy], xim (maneuvering) 
if m > 0 

% The function 'manmatrix' creates the matrix of maneuver times 
[mantime,xim] = manmatrix(A,B,t); 
x2m = zeros(size(xim)); 
for k = I :m 

ztrue = [ztrue;H*xim(:,k)]; % stacks the ztrues to make a 2*N x I vector 
x2m(:,k) = timestep(E,xim(:,k)); 

end 

% Update posout matrix with the stacked ztrues: 

%[xl;yl;x2;y2;[non-maneuvering x3;y3;x4;y4;...maneuvering] 
posout = [posout,ztrue]; 
end 

% INITIALIZATION: Assume all targets are moving in a straight hne! 

[rr,II] = size(mantime); 
simt = simt + delta; 

%Sp = H*x2m(:,I); % ASSUMES THAT THE SENSOR PLATEORM IS THE HRST 
ELEMENT OE M-SET 

aUxI = [xin, xim]; 
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allx2 = [x2n,x2m]; 

% First, grab the sensor information to get the apphcable errors: 

% Perform a Kalman filter initiahzation for every track - assume straight hne! 
ztemp = []; 
for k = 1 :ns 
ser = nsens(k); 

[Q,P,R,D,derror,xhat,zret] = init(aUxl,albc2,delta,Spi,Sp,ser); 
xhatl(:,:,k) = xhat; 

Rl(:,:,k) = R; 
zl(:,:,k) = zret; 

% For multiple sensors, the z's are stacked, not 2-D to begin with 

% [(tgtl,sensl);(tgt2,sensl);.(tgtl,sens2);(tgt2,sens2);...] 

Ql(:,:,k) = Q; 

Pl(:,:,k) = P; 

derrorl(:,:,k) = derror; % Same for the errors, one number per point 

% [(tgtl,sensl);(tgt2,sensl);.(tgtl,sens2);(tgt2,sens2);...] 

end 

% Call init for ownship: 

[spP,spR,sper,Spst,spmeas] = spinit(spxi,spx2,Q,D); 

% Create the noisy measurements matrix for output: 
spzout = [spzout,spmeas]; 

zout = [zout,zl]; % Noisy z measurements zl;z2., and ns deep (1 page for each sensor) 

zl = []; 

state = [state, xhat 1]; 
spstate = [spstate,Spst]; 

Sptrue = H*spx2; % true x-y for ownship 

Sp = l-I*Spst; % STATE predicted x-y for ownship 

error = [error,derrorl]; 

sperror = [sperror,sper]; 

derror 1 = []; 

% Set X equal to the latest position, 
xin = x2n; 
xim = x2m; 

ztrue = setztrue(m,n,xin,xim); 
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posout = [posout,ztrue]; 

% indx will be used to keep track of where the targets are in their maneuver cycle: 
indx = ones(l,m); 
in= 1; 

for ii = 3:nsamples 
% Iterate the targets: 

% NON-MANEUVERING TARGETS, AI, xin, 
if n > 0 
fork= I:n 

xin(:,k) = timestep(E,xin(:,k)); 
end 
end 

% MANEUVERING TARGETS, A, xim, xm 
if m > 0 
index = I; 

for k = I :m % must loop thru for each target before moving on 
% Check for zeros here? 
if indx(k) < rr % Executes every time 
LEG = mantime(indx(k),k); 

TURN = mantime(indx(k)-i-I,k); 
end 

if indx(k) == rr 
if mod(rr,2)==0 

LEG = mantime(indx(k),k); 

TURN = mantime(indx(k)-i-I,k); 
elseif mod(rr,2) == I 

LEG = mantime(indx(k),k); 

TURN = LEGh-I; 
end 
end 

if ((simt <= LEG) I (simt > TURN) I (TURN == 0)) 
xim(:,k) = timestep(E,xim(:,k)); % Use non-maneuvering E if tgt is SLM 
elseif (simt > LEG) & (simt <= TURN) I (LEG == 0) % What to do if turn is a zero? or 


Leg? 
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xim(:,k) = timestep(Fturn(index:(index+3),:),xim(:,k)); 
end 

% Reset LEG & TURN as needed for each target: 
if simt > TURN & TURN > 0 & indx(k) < rr 
indx(k) = indx(k) + 2; 
end 

index = index+4; 
end 
end 

measurement h^******************^^*** 

ztrue = setztrue(m,n,xin,xim); 
posout = [posout,ztrue]; 
zre= []; 

zre = reshape(ztrue,2,mn); 

% SENSOR PLATEORM MEASUREMENT AND TRACKING (NAV solution): 
sensor =12; 
if flag == 0 

spxtrue = timestep(E,spx2); 

[spR,spmeas,sper] = meas(sensor,spxtrue,[0;0]); 
elseif flag == 1 

[ spxtrue, spmeas, sper, spR,in] = spitmeas(sensor, spx2, spman,ownfturn, simt, in); 
end 

[Spst,count] = kalmanl(Q,spP,spR,Spst,spmeas,count); % Calc sensor position state 

spx2 = spxtrue; 

sppos = [sppos,H*spxtrue]; 

spstate = [spstate,Spst]; 

sperror = [sperror,sper]; 

spzout = [spzout, spmeas]; 

Sp = H*Spst; % The x-y position of the sensor platform (nav solution) 

% TARGET MEASUREMENT (for each sensor) 
for k = l:ns 
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sensor = nsens(k); % Sends the sensor type (numerical) to the measurement function 
zret = []; 

[R,zret,derror] = meas(sensor,zre,Sp); 

Rl(:,:,k)=R; 
zl(:,:,k) = zret; 
derrorl(:,:,k) = derror; 
end 

error = [error,derrorl]; 

% Create the noisy measurements matrix for output: 
zout = [zout,zl]; 

% zret are the measurement inputs for the tracker 
% Xhat is the matrix of updated states 

^ H^************^^***** tracker 

for k = l:ns 

[xhat33,count] = tracker(Ql(:,:,k),Pl(:,:,k),Rl(:,:,k),xhatl(:,:,k),zl(:,:,k),count,trkr(k)); 
xhatl(:,:,k) = xhat33; 
end 

state = [state,xhatl]; 

^ DATA LINK FUNCTION *********>]<*>]<**>]<*>]<*>]<>]<*>]<*>]<* 

% The current sensor state for surveillance radar is used as a starting point, we're 
% getting farther away from the truth with hnk reports. Simulate by adding more error 
% Also, this function is only called every 5 loops (3 seconds right now). 


if mod(ii,5) == 0 

[hnkvector,lstate] = hnkreport(xhatl(:,:,l),simt); 

% Istate is simply the state vector for the data hnk report 
% hnkstate = [hnkstate,lstate]; 

% Irep is in the form 

[x;y;vel;IFFl;IFF2;IFF3;IFF4;IFFC;emitters;time;reportingunit,x2;.] 

hnkout = [linkout,hnkvector]; 
end 

count = count + 1; 
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simt = simt + delta; 

% SIMPLOT function for plotting on the fly 
% could use this for plotting the fused result 
%simplot(mn,ns,xhatl,Spst); 
end 

% **** CALCULATE THE LEAST SQUARE ERROR: 

for k = 1 :ns 

LSerr = lsecalc(mn,state(:,:,k),posout,nsamples); 

LSerror(:,:,k) = LSerr; 
end 


ifkk== 1, 
merror = error; 
zoutmean = zout; 
statemean = state; 
spstatemean = spstate; 
sperrormean = sperror; 
spmeasmean = spzout; 
hnkmean = linkout; 

LSmean = LSerror; 
else 

merror = merror + error; 
zoutmean = zoutmean + zout; 
statemean = statemean + state; 
spstatemean = spstatemean + spstate; 
sperrormean = sperrormean + sperror; 
spmeasmean = spmeasmean + spzout; 
hnkmean = hnkmean + hnkout; 
LSmean = LSmean + LSerror; 
end 


merror = merror/nloops; 
zoutmean = zoutmean/nloops; 
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statemean = statemean/nloops; 
spstatemean = spstatemean/nloops; 
sperrormean = sperrormean/nloops; 
spmeasmean = spmeasmean/nloops; 
linkmean = linkmean/nloops; 

LSmean = LSmean/nloops; 

% 

hand = trackplot(posout,zoutmean,statemean,sppos,spstatemean,linkmean); 
handl = errorplots(delta,merror,LSmean, trkr); 


NAME 

start 

INPUT 

a dummy variable 

OUTPUT 

A,Al,delta, t,nloops,nsamples,Spi,flag,ownsens 

PURPOSE 

The start function reads the start.mat file to collect timing information about the 
simulation, initiates the matrices of targets, and initiates other needed matrices (F,H) 


function [A,Al,delta,t,nloops,nsamples,Spi,flag,ownsens] = start(dumb) 


% OUTPUT VARIABLES: 

% A Maneuvering target matrix 

% A1 Non-maneuvering target matrix 

% delta delta time 

% t total simulation time in seconds 

% nloops number of simulation runs (for error analysis) 

% nsamples number of time steps in the simulation run 
% Spi Initial position of the sensor platform 

% flag Indicator of whether or not the sensor platform is maneuvering 
% ownsens Vector of ownship sensor types 
global H F; 

% START READS THE SIMULATION RUN PARAMETERS 
simstart = textread('start.mat'); 

% simstart returns in the following fashion: [t(simtime), nloops, deltatime, DLRP, 
% placeholder, tracker selection 1 -4] 
t = simstart(l)*60; % To put time into seconds for the simulation 
nloops = simstart(2); 

% 

%t = simtime(dumb); 
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%t = t*60; % to put time in seconds 

% determine the time increment and the total number of samples 
uprate = simstart(3); 

[delta,nsamples] = deltatime(t,uprate); 

% Start by getting file data on aU targets: A1 returns as a matrix of targets - N x 17 
A1 = textread('tardat.mat'); 

A=[]; 

Spi = [A1(1,2);A1(1,4)]; % First target of A1 is always ownship!! (the sensor platform) 
ownsens = [Al(l,10) Al(l,ll) Al(l,12) Al(l,13)]; 

% Now, determine if the target is maneuvering or not by checking item 17 
% and creating a new matrix, A and deleting the maneuvering tgt rows from Al: 

[d,e] = size(Al); 

Atemp = []; 
for i = 1 :d 

if Al(i,17)>0 
A=[A;Al(i,:)]; 
ifi== 1 
flag = 1; 
end 

elseif Al(i,17) == 0 

Atemp = [Atemp;Al(i,:)]; 
end 
end 

Al = []; 

Al = Atemp; 

% Now A contains the targets that are maneuvering only, for use by the 'fturn' function 
% Stack the splat data on A. Sensor platform is always the first row (or column in fturn) 
% If the sensor platform is not maneuvering, the splat info gets added to Al 
H= [1,0,0,0; 

0 , 0 , 1 , 01 ; 

F=[l, delta, 0, 0; 

0, 1, 0, 0; 

0, 0, 1, delta; 

0, 0, 0, 1]; 
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NAME 

fturn 

INPUT 

A, delta 

OUTPUT 

Fturn, B 

PURPOSE 

The fturn function reads the start.mat file to collect timing information about the 
simulation, initiates the matrices of targets, and initiates other needed matrices (F,H) 


function [Fturn,B] = fturn(A,delta) 


% This function wiU create a matrix containing submatrices of Fturn's for each target 
% platform, including the sensor platform, if maneuvering. 

B = []; 

Btemp = []; 

Fturn = []; 

% Initialize state vector of positions and velocities based on input values: 

[nn,b] = size(A); 
for jj = l:nn 
if A(ii,16) ==0 
turng =1.1; 
elseif A(jj,16) == 1 
turng = -1.1; 
end 

hdg = A(jij,5)*pi/180.0; 

speed = A(j[i,6)*6076.0/3600.0; % speed in feet/sec 
xi = [A(ii,2);speed*sin(hdg);A(jij,4);speed*cos(hdg)]; 

Btemp = [turng;hdg;speed;xi]; 
g = turng; 

% Turn omega calculation 
w = g*32.174/speed; 
r = w*delta; 
s = sin(r); 
c = cos(r); 

Ft = [l,(s/w),0,((l-c)/w); 

0,c,0,-s; 

0,((l-c)/w),l,(s/w); 

0,s,0,c]; 

Fturn = [Fturn;Ft]; 

B = [B,Btemp]; 
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end 


NAME 

setuptrackers 

INPUT 

ns 

OUTPUT 

trkr 

PURPOSE 

This function sets up the vector of desired trackers by reading the primary requested 
tracker and then assigning another tracker for the other sensors. 


function trkr = setuptrackers(ns) 


% 1 = PDAF 

% 2 = IMM 

% 3 = Kaknan filter 

% 4 = Const Gain Kalman 

% 5 = Kalman filter #2 (for navigation) 

trki- = []; 

simstart = textread('start.mat'); 

primarytracker = simstart(6); % Reads the desired tracker for the primary radar tracking 
for k = 1 :ns 

trkr= [trkr,primarytracker]; 
if primarytracker == 1 
temp = 2; 

elseif primarytracker == 2 
temp = 3; 

elseif primarytracker == 3 
temp = 4; 

elseif primarytracker == 4 
temp = 1; 
end 

primarytracker = temp; 
end 
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NAME 

manmatrix 

INPUT 

A,B,t 

OUTPUT 

mantime,xim 

PURPOSE 

This function sets up the Maneuver Time Matrix based on the size of the largest 
maneuver vector and calls the maneuvertime function with the needed info. Based on the 
target vector and simulation time, determine the times for maneuvers and creates a vector 
to hold the cumulative times. For consistency, xim is returned, similar to the nonman 
function 


function [mantime,xkn] = manmatrix(A,B,t) 
time = t; 

[N,b] = size(A); 
mantime = []; 
mt= []; 
forkk= 1;N 

legt = round(A(kk,15)*60); 
turnt = round(A(kk,17)*60); 
mt = maneuvertime(time,legt,turnt); 
ifkk== 1 
r = length(mt); 
longestr = r; 
mantime = zeros(r,N); 
mantime(:,l) = mt; 
end 

ifkk> 1 

s = length(mt); 
r = longestr; 
if s > r 

mantime = [mantime;zeros((s-r),N)]; % Adds s-r rows of zeros to mantime 
mantime(:,kk) = mt; 
longestr = s; 
elseif s < r 

mantime(:,kk) = [mt;zeros((r-s),l)]; 
elseif s == r 

mantime(:,kk) = mt; 
end 
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end 


end 

% Peel off the x's 
xim = B(4:7,:); 


NAME 

maneuvertime 

ENPUT 

t,legt,turnt 

OUTPUT 

Mt 

PURPOSE 

This function determines if the target is maneuvering based on the times given and uses 
legt and turnt to calculate the time at which the target turns and stops turn: t enters the 
function in seconds, returns mt which is a vector of turn times for a maneuvering target, 
or else it returns an empty vector for a target that wiU not maneuver. 


function [mt] = maneuvertkne(t,legt,turnt) 


mt=[]; 


if turnt == 0 

dispCAre you nuts? This target is not maneuvering. Program Terminated!'); 
return 
end 

totaltime = 0; 
step = 0; 
k=l; 

while (totaltime < t) 
if (mod(k,2) == 1) 
step = legt; 

elseif (mod(k,2) == 0) 
step = turnt; 
end 

totaltime = totaltime + step; 
if totaltime > t 
totaltime = t; 
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end 


mt = [mt;totaltime]; 
k = k + 1; 
end 

[m,o] = size(mt); 

if (mod(m,2) == 1) % if number of elements is odd 
mt = [mt;0]; 
end 

% Returns the vector with the times for maneuver in the following format: 
% [time at which turnl executes; 

% time to return to straight hue 
% time to start next turn 
% .] 


NAME 

nonman 

ENPUT 

A1 

OUTPUT 

xin 

PURPOSE 

This function will extract the position and velocity data form the A1 matrix of non¬ 
maneuvering targets. Of course, this also includes fixed targets. Initializes state vector 
of positions and velocities based on input values: 


function [xin] = nonman(Al) 

[N,b] = size(Al); 
xin = zeros(size(4,N)); 
forjij= 1:N 

hdg = Al(ii,5)*pi/180.0; 
if Al(ii,6) >0 

speed = Al(ii,6)*6076.0/3600.0; % speed in feet/sec 
elseif A1 (jj,6) == 0 
speed = 0; 
end 

xin(l:4,jij) = [Al(j[i,2);speed*sin(hdg);Al(j[i,4);speed*cos(hdg)]; 
end 
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NAME 

ownxi 

INPUT 

ownship 

OUTPUT 

Spxi 

PURPOSE 

This function extracts the position and velocity data from the ownship vector, initializes 
state vector of positions and velocities based on input values 


function [spxi] = ownxi(ownship) 
hdg = ownship(5)*pi/180.0; 

speed = ownship(6)*6076.0/3600.0; % speed in feet/sec 
spxi = [ownship(2);speed*sin(hdg);ownship(4);speed*cos(hdg)]; 
end 


NAME 

ftum 

INPUT 

aUx 1, albc2, delta, Spi, Sp, sensor 

OUTPUT 

Q,P,R,D,derror,xhat,zret 

PURPOSE 

This function provides the initial gain for use by any tracking filter, assumes straight hne 
motion, and provides noisy measurements via pole2cart. 


function [Q,P,R,D,derror,xhat,zret] = init(aLbci,aUx2,deita,Spi,Sp,sensor) 

% Input variables; 

% H extraction matrix 

% F Straight line motion matrix 

% allxl initial x positions 

% allx2 incremented x positions 

% delta time delta 

% Sensor (vector) with apphcable sensor for senserror 

% 

% Output variables: 

% Q q matrix 

% P initial covariance 

% derror error matrix 

% xhat initial X hat matrix [x;vx;y;vy,x;vx;y;vy,.] one xhat for each target. 

% zret xy measurement vector (noisy) 
global H F; 

% INITIALIZATION STEPS: 
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P = []; % Size of P will be 4*N x 4 (a 4x4 for each target) 

y=[]; 
y2 = []; 
xhat = []; 
a = delta^2/2; 
b = delta^3/3; 

% time for some special purpose matrices: 

Q = [b, a, 0, 0; 
a, delta, 0, 0; 

0, 0, b, a; 

0, 0, a, delta]; 

D = [1,0, 0, 0; 

1/delta, 0, -1/delta, 0; 

0, 1, 0, 0; 

0 1/delta, 0, -1/delta]; 

zret = []; 
ztemp = []; 
derror = []; 
derror2 = []; 

R = []; 

Rinit = []; 
sigr = 0.0; 
sigb = 0.0; 

Sv=[]; 

[sigr,sigb,Sv] = senserror(sensor); 

[aa,bb] = size(allxl); 

% Measurement 1: Add noise to initial positions: 
for h = 1 :bb 

% Take a noisy measurement for aU the targets (man and non-man) 
[zcart,Rl,disterror] = pole2cart(H*aUxl(:,h),Spi,sigr,sigb,Sv); 
zret = [zret;zcart]; 
y = [y,zcart]; 

derror = [derror;disterror]; 

Rinit = [Rinit;Rl]; 
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end 

% Measurement 2: 
for h = 1 ;bb 

% Take a noisy measurement! 

[zcart,R2,disterror] = pole2cart(H*albc2(:,h),Sp,sigr,sigb,Sv); 
ztemp = [ztemp;zcart]; 
y2 = [y2,zcart]; 

R= [R;R2]; 

derror2 = [derror2;disterror]; 
end 

zret = [zret,ztemp]; 
derror = [derror,derror2]; 
y=[y;y2]; 
ytemp = zeros(4,l); 

forh= l:bb 

% vector yflipped contains the z's stacked z(2)/z(l) 
ytemp = y(:,h); 

yflipped = [ytemp(3);ytemp(4);ytemp(l);ytemp(2)]; 
plantnoise = Rinit(2*h-l:2*h,:)+H*inv(F)*Q*(inv(F))'*H'; 
cov = [R(2*h-l:2*h,:) zeros(2,2); 

zeros(2,2) plantnoise]; 

PI = D*cov*D'; 

P= [P;P1]; 
xhatl = D* yflipped; 

xhat = [xhat;xhatl]; % xhat is returned as a stacked vector of [x;vx;y;vy. 

end 


154 




NAME 

polelcart 

INPUT 

ztrue,sp,sr,sb,Sv 

OUTPUT 

z,R,disterror 

PURPOSE 

This function processes the measurements from the sensors by converting to polar coords 
and adding the supphed error, then converting back to Cartesian. 


function [z,R,disterror] = pole2cart(ztrue,sp,sr,sb,Sv) 

% Input Arguments 
% ztrue: True target position 

% sp: True sensor position 
% sr: Range Standard deviation error 

% sb: Bearing standard deviation error 

% sv; sensor variance matrix (2x2) 

% Outputs: 

% z: Output measurement with error in cart coords 

% R: Covariance Matrix of the meas in cart coords 

% disterror: Distance error between meas and true posit 

% Convert current position to polar coordinates and add sensor noise 

% Generating a noisy measurement from the sensor: 

z = ztrue-sp; % position relative to the sensor 

r = sqrt(z(l)^2 + z(2)^2); % Range from sensor 

b = atan2(z(2),z(l)); % Bearing from sensor 

r = r+sr*randn; % Range + sensor noise 

b = b+sb*randn; % Bearing + sensor noise 

% Convert the measurement to Cartesian coordinates 
z = [r*cos(b);r*sin(b)]+sp; 

% Now make the measurement covariance in Cartesian coords: 
fx = [cos(b) -r*sin(b);sin(b) r*cos(b)]; 

R = fx*Sv*fx'; 

% Compute the measurement error in Cartesian coordinates 

ztilde = ztrue - z; 

disterror = sqrt(ztilde'*ztilde); 
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NAME 

spinit 

INPUT 

spxi,spx2,Q,D 

OUTPUT 

spP,spR,sperror,spstate,spmeas 

PURPOSE 

This function provides the initial gain and state estimation for ownship navigation, 
assumes straight hne motion, and provides noisy measurements via pole2cart. 


function [spP,spR,sperror,spstate,spmeas] = spinit(spxi,spx2,Q,D) 

% Input variables: 

% spxl initial X positions 

% spx2 incremented x positions 

% Q Q matrix 

% D D matrix 

% Output variables: 

% spP initial covariance 

% spR sensor error covariance 

% sperror error matrix 

% spstate initial xhat matrix [x;vx;y;vy,x;vx;y;vy,.] one xhat for each target. 

% spmeas noisy platform measurements 
global H F; 

% INITIALIZATION STEPS: 

spP=[]; 

spR=[]; 

y=[]; 
y2 = []; 
spstate = []; 
spQ = 1000*Q; 

sensorposition = [0;0]; % Assumes that the origin of the nav solution is ownship eg 
spmeas = []; 
sperror = []; 
sensor =12; 

[sigr,sigb,Sv] = senserror( sensor); 

%sigr = 20.0; % Feet 

%sigb = 0.001; % Radians from ownship 

%Sv = diag([sigr^2;sigb^2]); 

% Measurement 1: Add noise to initial position: 

% Take a noisy measurement 
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[zcart 1 ,spRl ,disterror] = pole2cart(H*spxi,sensorposition,sigr,sigb,Sv); 
spmeas = [spmeas,zcartl]; 
sperror = [sperror,disterror]; 

% Measurement 2: 

% Take a noisy measurement 

[zcart2, spR2,disterror] = pole2cart(H* spx2, sensorposition, sigr, sigb, S v); 
spmeas = [spmeas,zcart2]; 
sperror = [sperror,disterror]; 

% vector yflipped contains the z's stacked z(2)/z(l) 
yflipped = [zcart2(l);zcart2(2);zcartl(l);zcartl(2)]; 
plantnoise = spRl + H*inv(F)*spQ*(inv(F))'*H'; 
cov = [spR2 zeros(2,2); 

zeros(2,2) plantnoise]; 
spP = D*cov*D'; 
spstate = D*yf]ipped; 


NAME 

setztrue 

INPUT 

m,n,xin,xim 

OUTPUT 

ztrue 

PURPOSE 

This function simply extracts the true x,y data from the maneuvering/non-maneuvering 
vectors. 


function ztrue = setztrue(m,n,xin,xim) 
global H; 
ztrue = []; 

% NON-MANEUVERING 
if n > 0 
fork= l:n 

ztrue = [ztrue;H*xin(;,k)]; 
end 
end 

% MANEUVERING 
ifm>0 
fork= l:m 


ztrue = [ztrue;H*xim(;,k)]; 
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end 

end 


NAME 

senserror 

INPUT 

sensor 

OUTPUT 

sigr,sigb,Sv 

PURPOSE 

This function matches the type of sensor to the corresponding errors. For each of the 
sensors on the sensor platform, the range and bearing errors are returned along with the 
sensor variance matrix. 


function [sigr,sigb,Sv] = senserror(sensor) 


% NONE OF THE FORMAT TYPES ARE CORRECTED!! ALL SENSORS ARE IN X-Y 
COORDS 

% Input variables: 

% sensor # representing the sensor type 
% I Airborne search radar 
% 2 Shipboard search radar 

% ... 

% See the switch-case statement for the rest of the descriptions. 

% Output variables: 

% sigr sigma range 
% sigb sigma bearing 
% Sv sensor variance matrix 
switch sensor 
case I, 

% AIRBORNE SURVEILLANCE RADAR 2-D 
sigr = 100; 
sigb = .005; 

Sv = diag([sigr^2;sigb^2]); 
case 2, 

% SURFACE SURVEILLANCE RADAR 2-D 
sigr = 20; 
sigb = .0001; 

Sv = diag([sigr^2;sigb^2]); 
case 3, 


158 








% AIRBORNE FIGHTER RADAR 
sigr = 20; 
sigh = .001; 

Sv = diag([sigr^2;sigb^2]); 
case 4, 

% SURFACE FIRE CONTROL RADAR 
sigr =10; 
sigb = .00001; 

Sv = diag([sigr^2;sigb^2]); 
case 5, 

% NCTR 
sigr = 20; 
sigb = .0001; 

Sv = diag([sigr^2;sigb^2]); 
case 6, 

% FIFE 
sigr = 20; 
sigb = .0001; 

Sv = diag([sigr^2;sigb^2]); 
case 7, 

% IFF INTERROGATE 
sigr = 50; 
sigb = .001; 

Sv = diag([sigr^2;sigb^2]); 
case 8, 

% IFF REPLY ONLY (NOT USED YET) 
sigr = 20; 
sigb = .0001; 

Sv = diag([sigr^2;sigb^2]); 
case 9, 

% AIRBORNE ESM 
sigr = 500; 
sigb = .01; 

Sv = diag([sigr^2;sigb^2]); 
case 10, 
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% SURFACE ESM 
sigr = 300; 
sigh = .01; 

Sv = diag([sigr^2;sigb^2]); 
case 11, 

% IR DETECTOR 
sigr = 6000; 
sigb = .0001; 

Sv = diag([sigr^2;sigb^2]); 
case 12, 

% OWNSHIP NAVIGATION!! 
sigr = 10; 
sigb = 0.0001; 

Sv = diag([sigr^2;sigb^2]); 
case 13, 

% DATA LINK REPORT TRACKING ERRORS 
sigr = 200; 
sigb = .01; 

Sv = diag([sigr^2;sigb^2]); 
otherwise, 

dispCUnknown sensor') 
end 


NAME 

timestep 

INPUT 

P,x 

OUTPUT 

xnew 

PURPOSE 

The timestep function receives the time increment (E-matrbc) and returns the new state 
vector for a moving target or platform, can be E or Eturn. 


function [xnew] = timestep(E,x) 


xnew = E*x; 
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NAME 

spitmeas 

INPUT 

sensor,spx2,spman,ownft urn,simt,in 

OUTPUT 

spxtrue, spmeas, sper, spR,in 

PURPOSE 

This function will iterate the sensor platform and calculate the noisy measurement 


function [spxtrue,spmeas,sper,spR,in] = spitmeas(sensor,spx2,spman,ownfturn,smit,m); 


global H F sLEG sTURN; 

origin = [0;0]; % Assumes that the origin of the nav solution is ownship eg 

spmeas = []; 

sper= []; 

spR=[]; 

sLEG = 0; 

sTURN = 0; 

% spman has the maneuver times, only need to do one test: 

[rr,ll] = size(spman); % so rr = the number of maneuver times 
% Iterate the truth position: 

% Check for zeros here? 
if in < rr % Executes every time 
sLEG = spman(in); 
sTURN = spman(in+l); 
end 

if in == rr 

ifmod(rr,2)==0 
sLEG = spman(in); 
sTURN = spman(in+l); 
elseif mod(rr,2) == 1 
sLEG = spman(in); 
sTURN = sLEG+1; 
end 
end 

if ((simt <= sLEG) I (simt > sTURN) I (sTURN == 0)) 

spxtrue = timestep(E,spx2); % Use non-maneuvering E if tgt is SLM 
elseif (simt > sLEG) & (simt <= sTURN) I (sLEG == 0) % What to do if turn is a zero? 

or Leg? 

spxtrue = timestep(ownfturn,spx2); 
end 
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% Reset LEG & TURN as needed for each target; 
if simt > sTURN & sTURN > 0 & in < rr 
in = in + 2; 
end 

% Measurement: Add noise to position: 

[sigr,sigb,Sv] = senserror( sensor); 

[spmeas,spR,sper] = pole2cart(H*spx2,origin,sigr,sigh,Sv); 


NAME 

meas 

INPUT 

sensor,zre,Sp 

OUTPUT 

R,zret,derror 

PURPOSE 

This function takes the noisy measurements and returns the matrix of noisy 
measurements (depending on the sensor type) 


function [R,zret,derror] = meas(sensor,zre,Sp) 


% INPUTS: 

% sensor: sensor type (just search radar for now) 

% zre: x matrix of true x,y measurements 
% Sp: sensor platform position 
% OUTPUTS: 

% xret output matrix of noisy x,y coordinates 
% disterror actual distance errors for each measurement 

% R error covariance for the measurement of each target [R1;R2;R3;R4;R5....] 

%global MEAS_R zpol; 

zret= []; 

derror = []; 

sigr = 0.0; 

sigb = 0.0; 

Sv=[]; 

[sigr, sigb, Sv] = senserror( sensor); 

[aa,bb] = size(zre); 

R = zeros(2*bb,2); 
forh= l:bb 
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% Take a noisy measurement! 

[zcart,Rl,disterror] = pole2cart(zre(:,h),Sp,sigr,sigb,Sv); 
zret = [zret;zcart]; 
derror = [derror;disterror]; 

R(2*h-l:2*h,:) = Rl; 
end 


NAME 

tracker 

INPUT 

Ql,Pl,Rl,xhatl,zl,count, trkr 

OUTPUT 

xhat33,count 

PURPOSE 

Right now, this function only calls the appropriate tracker with aU the apphcable track 
data. 


function [xhat33,count] = tracker(Ql,Pl,Rl,xhatl,zl,count,trkr) 
switch trkr 
case 1, 

[xhat33,count] = pdaf(Ql,Pl,Rl,xhatl,zl,count); 

case 2, 

[xhat33,count] = immtrk(Ql,Pl,Rl,xhatl,zl,count); 
case 3, 

[xhat33,count] = kalman(Ql,PI,Rl,xhatl,zl,count); 
case 4, 

[xhat33,count] = constGKF(Ql,Pl,Rl,xhatl,zl,count); 
case 5, 

[xhat33,count] = kalmanl(Ql,PI,Rl,xhatl,zl,count); 
otherwise, 

dispCUndefined Tracker') 
end 
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NAME 

pdaf 

INPUT 

Q,P,R,xhat,zret,count 

OUTPUT 

xhat33,count 

PURPOSE 

Tracker for fusiin, will perform the following functions: 

PREDICT, based on previous stuff (stored as globals) 

UPDATE, based on current measurement 


% Thesis Tracker: Probabilistic Data Association Filter 
binction [xhat33,count] = pdaf(Q,P,R,xhat,zret,count) 
global H F M I Qp Pd Pg Pp lambda sig2r sig2b gatev; 
global Sp; 
if count > 1 

xhatone = zeros(size(xhat)); 
m = 2; 
fork= 1:M 

a = 2*k-l; 
b = 2*k; 
c = 4*k-3; 
d = 4*k; 
xhatest = []; 

Pproj = []; 

% Now for the weighted predictions: 
xhatone(c:d,l) = F*xhat(c:d,l); 

Pproj = F*Pp(c:d,:)*F' + Qp(c:d,:); %okay so far 

% Calculate the b using non-parametric assumptions (diffuse prior model) 

e=[]; 

nu = []; 

nuk = zeros(2,l); 

beta = []; 

mess = zeros(2,2); 


% Semi-major and semi-minor axis in feet: 
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% minor = 600; 

%major = 1500; 

%Vk = 2*pi*minor*major; 

% lambda = m/Vk; 

PreS = H*Pproj*H'; 

% Calculate the innovation stuff for the actual measurement, based on meas S: 
Sk = PreS + R(a:b,:); 

Sinv = inv(Sk); 

K = Pproj*H'*Sinv; % W!! 

W = K; 

bb = lambda*((det(2*pi*Sk))^.5)*(l-Pd*Pg)/Pd; 


sum = bb; 

knu = zret(a:b,l) - H*xhatone(c:d,l); 

%knu = zcart - zproj; 

e(l) = exp(-0.5*knu'*Sinv*knu); 
nu(:,l) = knu; 
sum = sum+e(l); 

if m > 1, 
for j = 2:m 

[zclutter,Rclut,er] = pole2cart(zret(a:b, l),Sp,sig2r,sig2b,gatev); 
zset = zclutter; 

knu = zset - H*xhatone(c:d,l); 

e(j) = exp(-0.5*knu'*Sinv*knu); 
nu(:,j) = knu; 
sum = sum+e(j); 

end 

end 


for j = l:m 
v = nu(;,j); 
betafj) = e(j)/sum; 
nuk = nuk + betafj) *v; 
end 
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for j = l:m 
v = nu(:,j); 

mess = mess + (beta(j)*v*v'); 
end 

mess = mess-nuk*nuk'; 

betazero = bb/sum; 

Ptilde = W*mess*W'; 


Pc = Pproj-W*Sk*W'; 

Pp(c:d,:) = betazero*Pproj + (1-betazero)* Pc + Ptilde; 
xhatest = xhatone(c;d,l) + W*nuk; 
xhat33(c:d,l) = xhatest; 
end 
return 

elseif count == 1 

Pd = 1.0; 

Pg = .99; 
m = 2; 

% Used only for generating "noisier" points: 
sig2r = 900; 
sig2b = 2.0*pi/180; 
gatev = [sig2r^2,0; 

0,sig2b^2]; 

qsquared = 1000; 

I = eye(4); 

Pp = zeros(size(P)); 

Pp(:,:) = P; 

M = max(srze(xhat)); 

M = M/4; 

Qp = zeros(srze(P)); 
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xhatone = zeros(size(xhat)); 


for k = 1:M 

a = 2*k-l; 
b = 2*k; 
c = 4*k-3; 
d = 4*k; 

Qp(c:d,:) = qsquared*Q; 
xhatest = []; 

Pproj = []; 

% Now for the weighted predictions: 
xhatone(c:d,l) = F*xhat(c:d,l); 

Pproj = F*Pp(c:d,:)*F' + Qp(c:d,:); %okay so far 

% Calculate the b using non-parametric assumptions (diffuse prior model) 

e=[]; 

nu = []; 

nuk = zeros(2,l); 

beta = []; 

mess = zeros(2,2); 

% Semi-major and semi-minor axis in feet: 
minor = 600; 
major = 1500; 

Vk = 2*pi*minor*major; 
lambda = m/Vk; 

PreS = l-I*Pproj*l-l'; 

% Calculate the innovation stuff for the actual measurement, based on meas S: 
Sk = PreSH-R(a:b,:); 

Sinv = inv(Sk); 

K = Pproj*H'*Sinv; % W!! 

W = K; 

bb = lambda*((det(2*pi*Sk))^.5)*(l-Pd*Pg)/Pd; 
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sum = bb; 

knu = zret(a:b,l) - H*xhatone(c;d,l); 
e(l) = exp(-0.5*knu'*Sinv*knu); 
nu(:,l) = knu; 
sum = sum+e(l); 

if m > 1, 
for j = 2:m 

[zclutter,Rclut,er] = pole2cart(zret(a:b, l),Sp,sig2r,sig2b,gatev); 
zset = zclutter; 

knu = zset - H*xhatone(c:d,l); 

e(j) = exp(-0.5*knu'*Sinv*knu); 
nu(:,j) = knu; 
sum = sum+e(j); 

end 

end 


for j = l:m 
v = nu(;,j); 
beta(j) = e(j)/sum; 
nuk = nuk + beta(j)*v; 
end 

for j = l;m 
v = nu(;,j); 

mess = mess + (beta(j)*v*v'); 
end 

mess = mess-nuk*nuk'; 

betazero = bb/sum; 

Ptilde = W*mess*W'; 


Pc = Pproj-W*Sk*W'; 

Pp(c:d,:) = betazero*Pproj + (l-betazero)*Pc + Ptilde; 
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xhatest = xhatone(c:d,l) + W*nuk; 
xhat33(c:d,l) = xhatest; 

end 

return 


NAME 

immtrk 

INPUT 

Q,P,R,xhat,zret,count 

OUTPUT 

xhat33,count 

PURPOSE 

Tracker for fusim, will perform the following functions: 

PREDICT, based on previous stuff (stored as globals) 

UPDATE, based on current measurement 


% Thesis Tracker: Interacting Multiple Models (2) Filter 

function [xhat33,count] = knmtrk(Q,P,R,xhat,zret,count) 

global H F M I Qi Qturn Pi Pturn xhatone xhattwo alpha beta row musave; 


if count > 1 

for k = 1:M 

% PREDICTION STEP: 

a = 2*k-l; 
b = 2*k; 
c = 4*k-3; 
d = 4*k; 

mu = musave(:,k); 

% Pre-process the cbars and mu's: 

cbarone = row(l,l)*mu(l)H-row(2,l)*mu(2); 

cbartwo = row(l,2)*mu(l)-i-row(2,2)*mu(2); 

% Set up the calculations for xhat-zero-one 

xhatzone = xhatone(c:d,l)*(row(l,l)*mu(l)/cbarone) -H 
xhattwo(c:d, l)*(row(2, l)*mu(2)/cbarone); 

xhatztwo = xhatone(c:d,l)*(row(l,2)*mu(l)/cbartwo) + 
xhattwo(c: d, 1) * (row(2,2) *mu(2)/cbartwo); 

muoneone = row(l,l)*mu(l)/cbarone; 
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mutwoone = row(2,l)*mu(2)/cbarone; 
muonetwo = row(l,2)*mu(l)/cbartwo; 
mutwotwo = row(2,2)*mu(2)/cbartwo; 

xtildell = xhatone(c:d,l)-xhatzone; 
xtilde21 = xhattwo(c:d,l)-xhatzone; 
xtildel2 = xhatone(c:d,l)-xhatztwo; 
xtilde22 = xhattwo(c:d,l)-xhatztwo; 

Pzeroone = muoneone*(Pi(c:d,:)+xtildel l*xtildel 1') + 
mutwoone*(Pturn(c:d, :)+xtilde21 *xtilde21'); 

Pzerotwo = muonetwo*(Pi(c:d,;)+xtildel2*xtildel2') + 
mutwotwo * (Pturn(c; d,:) +xtilde22 * xtilde22'); 

% Now for the weighted predictions: 
xhatone(c:d,l) = F*xhatzone; 
xhattwo(c:d,l) = F*xhatztwo; 

Pone = F*Pzeroone*F' + Qi(c:d,:); 

Ptwo = F*Pzerotwo*F' + Qturn(c:d,:); 

% Kalman Straight hne: 

K = Pone*l-I'*inv(H*Pone*FI'+R(a:b,:)); 

% Kalman Turn: 

Kturn = Ptwo*l-l'*inv(H*Ptwo*l-I'+R(a:b,:)); 

Pi(c:d,:) = (I-K*H)*Pone*(I-K*H)'+K*R(a:b,:)*K'; 

Pturn(c:d,:) = (I-Kturn*H)*Ptwo*(I-Kturn*H)'+Kturn*R(a:b,:)*Kturn'; 

51 = H*Pone*H’+R(a:b,:); 

52 = H*Ptwo*H’+ R(a:b,:); 

% MEASUREMENT UPDATE 

% Calculate the state estimate for straight hne and turning: 

% Straight line 

ztl =zret(a:b,l) - l-l*xhatone(c:d,l); 
xhatest = xhatone(c:d,l)+K*(ztl); 

% Turning 

zt2 = zret(a:b,l) - H*xhattwo(c:d,l); 
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xhatturn = xhattwo(c;d,l)+Kturn*(zt2); 


% Scores: 

%m = 2.0; 

Lambdal = exp(-((ztr*inv(Sl)*ztl)/2.0))/( 2*pi*norm(S 1)^.5 ); 

Lambda2 = exp(-((zt2'*inv(S2)*zt2)/2.0))/( 2*pi*norm(S2)^.5 ); 
see = Lambdal *cbarone + Lambda2*cbartwo; 

%Store for next iteration 
xhatone(c:d,l) = xhatest; 
xhattwo(c:d,l) = xhatturn; 

% Combine the data for output 

xhat33(c:d,l) = mu(l)*xhatest + mu(2)*xhatturn; 

mustr = Lambdal *cbarone/see; 
muturn = Lambda2*cbartwo/see; 
musave(:,k) = [mustr;muturn]; 
end 

%musave 
elseif count == 1 

% Set the probability of a turn if in SLM (alpha) and the prob of stop turn (beta) 

% if in a turn: 

alpha = 0.1; 

beta = 0.33333; 

qsquared = 1; 

turnq2 = 100000; 

row = [(1-alpha),alpha;beta,(1-beta)]; 
xhatone = xhat; 
xhattwo = xhatone; 

1 = eye(4); 

Pi = zeros(size(P)); 

Pi(:,:) = P; 

Pturn = zeros(size(P)); 

Pturn(:,:) = P; 

M = max(size(xhat)); 
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M = M/4; 

% Initial State Likelihood: (straight hne track) 
musave = zeros) size(2,M)); 
mu = [1;0]; 

Qi = zeros(size(P)); 

Qturn = zeros(srze(P)); 
fork= 1:M 

musave(l:2,k) = mu; 
end 

fork= 1:M 

% INITIAL PREDICTION STEP ONLY!!: 


a = 2*k-l; 
b = 2*k; 
c = 4*k-3; 
d = 4*k; 

Qi(c:d,:) = qsquared*Q; 

Qturn(c:d,:) = turnq2*Q; 

% Pre-process the cbars and mu's: 

cbarone = row(l,l)*mu(l)-i-row(2,l)*mu(2); 

cbartwo = row(l,2)*mu(l)H-row(2,2)*mu(2); 

% Set up the calculations for xhat-zero-one 

xhatzone = xhatone(c:d,l)*(row(l,l)*mu(l)/cbarone) H- 
xhattwo(c: d, 1) * (row(2,1) *mu(2)/cbarone); 

xhatztwo = xhatone(c:d,l)*(row(l,2)*mu(l)/cbartwo) H- 
xhattwo(c: d, 1) * (row(2,2) *mu(2)/cbartwo); 

muoneone = row(l,l)*mu(l)/cbarone; 
mutwoone = row(2,l)*mu(2)/cbarone; 
muonetwo = row(l,2)*mu(l)/cbartwo; 
mutwotwo = row(2,2)*mu(2)/cbartwo; 

xtildell = xhatone(c:d,l)-xhatzone; 
xtilde21 = xhattwo(c:d,l)-xhatzone; 
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xtildel2 = xhatone(c;d,l)-xhatztwo; 
xtilde22 = xhattwo(c;d,l)-xhatztwo; 

Pzeroone = muoneone*(Pi(c;d,:)+xtildel l*xtildel 1') + 
mutwoone*(Pturn(c:d, :)+xtilde21 *xtilde21'); 

Pzerotwo = muonetwo*(Pi(c:d,:)+xtildel2*xtildel2') + 
mutwotwo * (Pturn(c; d,:) +xtilde22 * xtilde22'); 

% Now for the weighted predictions: 
xhatone(c:d,l) = F*xhatzone; 
xhattwo(c:d,l) = F*xhatztwo; 

Pone = F*Pzeroone*F' + Qi(c:d,:); 

Ptwo = F*Pzerotwo*F' + Qturn(c:d,:); 

% Kalman Straight hne: 

K = Pone*H'*inv(H*Pone*H'+R(a:b,:)); 

% Kalman Turn: 

Kturn = Ptwo*H'*inv(H*Ptwo*H'+R(a:b,:)); 

Pi(c:d,:) = (I-K*H)*Pone*(I-K*H)'+K*R(a:b,:)*K'; 

Pturn(c:d,:) = (I-Kturn*H)*Ptwo*(I-Kturn*H)'+Kturn*R(a:b,:)*Kturn'; 

51 =H*Pone*H’+ R(a:b,:); 

52 = H*Ptwo*H'+ R(a:b,:); 

% MEASUREMENT UPDATE 

% Calculate the state estimate for straight hne and turning: 

% Straight hne 

ztl = zret(a:b,l) - H*xhatone(c:d,l); 
xhatest = xhatone(c:d,l)+K*(ztl); 

% Turning 

zt2 = zret(a:b,l) - l-t*xhattwo(c:d,l); 
xhatturn = xhattwo(c:d,l)+Kturn*(zt2); 

xhat33(c:d,l) = mu(l)*xhatest + mu(2)*xhatturn; 

% Scores: 

%m = 2.0; 
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Lambdal = exp(-((ztr*inv(Sl)*ztl)/2.0))/( 2*pi*norm(S 1)^.5 ); %(2*pi)^(m/2)*norm(S 1)^.5 

Lambda2 = exp(-((zt2'*inv(S2)*zt2)/2.0))/( 2*pi*norm(S2)^.5 ); %(2*pi)^(m/2)*norm(S2)^.5 

see = Lambdal *cbarone + Lambda2*cbartwo; 

mustr = Lambdal *cbarone/see; 

muturn = Lambda2*cbartwo/see; 

musave(;,k) = [mustr;muturn]; 

xhatone(c:d,l) = xhatest; 
xhattwo(c:d,l) = xhatturn; 

end 

return 


NAME 

kalman 

INPUT 

Q,P,R,xhat,zret,count 

OUTPUT 

xhat33,count 

PURPOSE 

Tracker for fusim, will perform the following functions: 

PREDICT, based on previous stuff (stored as globals) 

UPDATE, based on current measurement 


% Thesis Tracker: Kalman Filter 

function [xhat33,count] = kalman(Q,P,R,xhat,zret,count) 

global H F Mk 1 Pk Qk; 


if count > 1 

xhat32 = zeros(size(xhat)); 
fork= LMk 

% PREDICTION STEP: 
Ptemp= []; 
a = 2*k-l; 
b = 2*k; 
c = 4*k-3; 
d = 4*k; 

Ptemp = (F*Pk(c:d,:)*F' + Qk(c:d,:)); 
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Pk(c:d,:) = Ptemp; 
xhat32(c:d,l) = F*xhat(c:d,l); 

% UPDATE STEP; 

K = Pk(c:d,:)*H'*inv(H*Pk(c:d,:)*H' + R(a:b,:)); % Slv = R 
Ppre(c:d,:) = (I-K*H)*Pk(c;d,:)*(I-K*H)'+K*R(a:b,:)*K'; 
xhat33(c;d,l) = xhat32(c;d,l)+K*(zret(a:b,l) - H*xhat32(c:d,l)); 
end 

Pk = Ppre; 

elseif count == 1 
qsquared = 10000; 

K=[]; 

I = eye(4); 

xhat32 = zeros(size(xhat)); 

Pk = zeros(size(P)); 

Mk = max(size(xhat)); 

Mk = Mk/4; 

Qk = zeros(size(P)); 
for k = l;Mk 

% INITIAL PREDICTION STEP ONLY!!: 

a = 2*k-l; 
b = 2*k; 
c = 4*k-3; 
d = 4*k; 

Qk(c:d,:) = qsquared*Q; 

Pk(c;d,:) = E*P(c:d,;)*E' + Qk(c;d,:); 
xhat32(c;d,l) = E*xhat(c:d,l); 

% UPDATE STEP: 

K = Pk(c:d,:)*H'*inv(H*Pk(c:d,:)*H' + R(a:b,:)); % Slv = R 
Ppre(c;d,:) = (I-K*H)*Pk(c;d,:)*(I-K*H)’+K*R(a:b,:)*K'; 
xhat33(c:d,l) = xhat32(c:d,l)+K*(zret(a;b,l) - H*xhat32(c:d,l)); 
end 
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Pk = Ppre; 


return 

end 


NAME 

constGKF 

INPUT 

Q,P,R,xhat,zret,count 

OUTPUT 

xhat33,count 

PURPOSE 

Tracker for fusim, will perform the following functions: 

PREDICT, based on previous stuff (stored as globals) 

UPDATE, based on current measurement 


% CONSTANT GAIN KALMAN FILTER 
function [xhat33,count] = constGKF(Q,P,R,xhat,zret,count) 
global Kbar M F H; 
if count > 1 

xhat32 = zeros(size(xhat)); 
for k = LM 

% PREDICTION STEP: 

a = 2*k-l; 
b = 2*k; 
c = 4*k-3; 
d = 4*k; 

xhat32(c:d,l) = F*xhat(c:d,l); 

% UPDATE STEP: 

xhat33(c:d,l) = xhat32(c:d,l)+Kbar(c:d,:)*( zret(a:b,l) - H*xhat32(c:d,l)); 
end 

elseif count == 1 
r = zeros(2,2); 

Kbar= []; 

qsquared = 10000.0; 

Q = qsquared*Q; 
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xhat32 = zeros(size(xhat)); 

M = max(size(xhat)); 

M = M/4; 
fork= 1:M 

% PREDICTION STEP: 

a = 2*k-I; 
b = 2*k; 
c = 4*k-3; 
d = 4*k; 

r(I:2,I:2) = R(a:b,I:2); 

G = eye(4); 

[KBAR,Pbar] = dIqe(F,G,H,Q,r); 

Kbar = [Kbar;KBAR]; 
xhat32(c:d,I) = F*xhat(c:d,I); 

% UPDATE STEP: 

xhat33(c:d,I) = xhat32(c:d,I)+KBAR*(zret(a:b,I) - H*xhat32(c:d,I)); 
end 
end 


NAME 

linkreport 

EVPUT 

xhat,simt 

OUTPUT 

linkdat,Istate 

PURPOSE 

This function creates a simulated data link report for the targets and allows for addition 
of attribute data, such as IFF data FORMAT: 

[x, y,velocity,heading,time,rffI,iff2,ifG,rff4,iffc,emitter I,emitter2,ID,reportingunit]; 


function [knkdatjstate] = Iinkreport(xhat,smit) 
global M; 
linkdat = []; 

Istate = []; 

%tempxhat = [0;0]; 
fork= I:M 
a = 4*k-3; %I 
b = 4*k-2; %2 
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c = 4*k-l; %3 
d = 4*k; %4 
% tempxhat(l) = xhat(a); 

% tempxhat(2) = xhat(c); 

% Add in noise 
X = xhat(a) + randn*50; 
y = xhat(c) + randn*50; 
vx = xhat(b) + randn*10; 
vy = xhat(d) + randn*10; 

% Add in a set value of bias error, simulates DLRP/gridlock errors 
X = X + 100; % 100 feet in bias error x-direction 
y = y + 100; % 100 feet in bias error y-direction 
hdg = 0; %round(atan(vx,vy)*180/pi); 

V = sqrt(vx^2 + vy^2); 

linkvector = [x, vx,y, vy, v,hdg,simt,0,0,0,0,0,0]; 

Istate = [lstate;linkvector]; 
hnkdat = [linkdat;x;y]; 
end 


NAME 

simplot 

INPUT 

mn,ns,xhat,Spst 

OUTPUT 

Plot with coordinates in feet 

PURPOSE 

This function is used to plot the main display (for each of the targets) Ownship is the 
first target. 


function simplot(mn,ns,xhat,Spst) 

% mn - Total number of targets 
% xhat - Current tracking state, 
global H; 

xyplot = zeros(2,l); 
figure(l) 

axis([-250000 850000 -250000 1000000]); 

xlabel('LONGITUDE - in feet for now'); ylabel('LATITUDE - in feet for now'); 
title)['Single Sensor Tracking Solutions']); 
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hold on; 
pause(.Ol) 
handle2 = []; 
r=0; g=l;b=0; 

handle 1 = plot( Spst( 1), Spst(3), 'm*'erasemode', 'xor', 'markersize',6); 

% set(handle 1,'color', [r g b]) 

for k = 1 ;ns 
ifk== 1 
r=0; g=0; b=l; 
elseif k == 2 
r=l;g=0;b=0; 
elseif k == 3 
r=.75; g=.5; b=0; 
elseif k == 4 
r=0; g=.5; b=.5; 
end 

for hh = l:mn 

xyplot = H*xhat(4*hh-3:4*hh,:,k); 

handle2(hh) = plot(xyplot( 1) ,xyplot(2),'+', 'erasemode', 'xor', 'markersize',6); 
set(handle2(hh),'color',[r g b]) 

%plot(xyplot(l),xyplot(2),'*'); 

%set(handle2,'xposition',xyplot(l),'yposition',xyplot(2)); 

end 

end 
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NAME 

Isecalc 

INPUT 

mn,state,posout,nsamples 

OUTPUT 

LSE 

PURPOSE 

This function calculates the Least Square Error comparing the state vector to the truth 
position. 


function LSE = lsecalc(nin,state,posout,nsamples) 


LSE = []; 
for hh = l:mn 
if =4*hh-3; 
i2 = 4*hh-l; 
tempi = [state(il,:); 
state(i2,:)]; 

tempposout = [posout(2*hh-l,2:(nsamples));posout(2*hh,2;(nsamples))]; 
Kaknanerror = tempposout-templ; 
temp2= []; 
temp3 = 0.0; 

LSerror = []; 

for g = l:nsamples-l, 

temp2 = [Kalmanerror(l,g);Kalmanerror(2,g)]; 
temp3 = sqrt( temp2(l)^2 + temp2(2)^2 ); 

LSerror = [LSerror, temp3]; 

end 

LSE = [LSE;LSerror]; 
end 


180 








NAME 

errorplots 

INPUT 

delta,merror,LSmean,trkr 

OUTPUT 

handl, error plots for each target and/or sensor 

PURPOSE 

This function plots the mean measurement errors and least squares tracking errors for 
each sensor and target 


function handl = errorplots(delta,merror,LSmean,trkr) 

% merror is mean measurement error 
% statemean is the matrix of mean states 
% posout is the truth data 
% delta is the time step 
% nothing returned except the plot handle. 

[a b c] = size(LSmean); 
nsamples = max(size(merror)); 
time = [0;max(size(LSmean)-l)]*delta; 
merror = merror(;,2;end,:); 

for hh = l:a % Plot for each target: 
figure 

title(['Mean Error (-) vs. Time, 1 Sensor Type/4 Trackers/100 run(s)']); 
axis([0 900-1 1]); 

xlabel('Time, (seconds)'); ylabel('Error, (feet)'); 
hold on; 

%for k = 1 :c 

plot(time,LSmean(hh,:,l),'b-'); % PDA 
plot(time,LSmean(hh,:,2),'g-'); % IMM 
plot(time,LSmean(hh,:,3),'r-'); % KALMAN 
plot(time,LSmean(hh,:,4),'m-'); % CGK 
plot(time,merror(hh,:,l),'c-'); % Mean Error 
%end 

end 


%for k = 1 :c % Loop for each sensor 
% merror2 = merror(:,2:nsamples,k); 

% trkrID = trkr(k); 

% forhh=l:a % Loop for each target 
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% figure 

% plot(tmie,merror2(hh,;),'-',tmie,LSmean(hh,;,k),'-'); 

% iftrkrID==l 


Error (-)/Mean PDAF Error (-) vs. Time (1 run(s))']); 

Error (-)/Mean IMM Error (-) vs. Time (1 run(s))']); 

Error (-)/Mean Kalman Error (-) vs. Time (1 run(s))']); 

Error (-)/Mean Const Gain Kalman Error (-) vs. Time (1 run(s))']); 

% end 
% end 
%end 
handl = 0; 


% title(['Mean Dist 

% elseif trkrID == 2 

% title( ['Mean Dist 

% elseif trkrID == 3 

% title( ['Mean Dist 

% elseif trkrID == 4 

% title( ['Mean Dist 


NAME 

trackplot 

INPUT 

posout,zoutmean,statemean,sppos,spstatemean,linkmean 

OUTPUT 

hand, track plot 

PURPOSE 

This function plots the true position, the noisy measurements and the track states for each 
track and each sensor, and provides a single plot for comparison of aU sensors and 
multiple plots of the primary tracker 


function hand = trackplot(posout,zoutmean,statemean,sppos,spstatemean,linkmean) 
[mn lots ns] = size(statemean); 
mn = mn/4; 

%figure 
%hold on 

%pol = zeros(2,max(size(posout))); ******* 

— po 1" 

%for i = l:mn 

% pol(:,:) = posout(2*i-l:2*i,:); 

% plot(pol(l,;),pol(2, 

%end 

%title(['True Target Motion/Mean Measured Position (1 run(s))']); 
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%for j = l:ns 
% for i = l;mn 

% po2(:,:) = zoutmean(2*i-l:2*i,:,j); 

% plot(po2(l,:),po2(2, 

% end 
%end 

figure 
hold on 
Hline= []; 
hnecolor = [001]; 
for j = 1 :ns % For each sensor; 
if j == 1 
r=0; g=0; b=l; 
elseif j == 2 
r=l;g=0;b=0; 
elseif j == 3 
r=.75; g=.5; b=0; 
elseif j == 4 
r=0; g=.5; b=.5; 
end 

for i = 1 ;mn % Loop for each target: 
poo = statemean(4*i-3,:,j); % X values 
pol = statemean(4*i-l,:,j); % Y values 
po2 = [poo;pol]; 

Hhne(i) = plot(po2(l,:),po2(2, 

end 

set(Hhne, 'color', linecolor); 
hnecolor = [r g b]; 
end 

title(['Mean State (Tracking Solution) for AU Sensors(100 runs)']); 
%figure 

axis([-250000 850000 -250000 1000000]); 
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xlabel('LONGITUDE (feet)'); ylabel('LATITUDE (feet)'); 

%hold on 

plot(sppos(l,:),sppos(2,:),'-',spstatemean(l,:),spstatemean(3,:),'-'); 
%poc = zeros(2,max(size(statemean))); 
pol = zeros(2,max(size(posout)-l)); 
for i = l:mn % Plot actual for each target 
pol(:,:) = posout(2*i-l:2*i,2:end); 

% poc(:,:) = [statemean(4*i-3,:,l);statemean(4*i-l,:,l)]; 

plot(pol(l,:),pol(2,:),'-'); %,poc(l,:),poc(2,:),'-'); 
end 

% Plot the data link result: 
for i = l:mn 

link(:,:) = hnkmean(2*i-l:2*i,l:end); 
plot(hnk( 1,: ),hnk(2,:),'m-'); 
end 

%title(['Mean Tracker State/Link Report vs. True Position(l run)']); 
hand = 0; 
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